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1.0 INTRODUCTION 


1.1 Background 

The EOS ICESat mission is scheduled for launch on July 2001. Three 
major science objectives of this mission are: (1) to measure long-term changes in the 
volumes (and mass) of the Greenland and Antarctic ice sheets with sufficient 
accuracy to assess their impact on global sea level, and to measure seasonal and 
interannual variability of the surface elevation, (2) to make topographic 
measurements of the Earth's land surface to provide ground control points for 
topographic maps and digital elevation models, and to detect topographic change, and 
(3) to measure the vertical structure and magnitude of cloud and aerosol parameters 
that are important for the radiative balance of the Earth-atmosphere system, and 
directly measure the height of atmospheric transition layers. The spacecraft features 
the Geoscience Laser Altimeter System (GLAS), which will measure a laser pulse 
round-trip time of flight, emitted by the spacecraft and reflected by the ice sheet or 
land surface. This laser altimeter measurement provides height of the GLAS 
instrument above the ice sheet. The geocentric height of the ice surface is computed 
by differencing the altimeter measurement from the satellite height, which is 
computed from Precision Orbit Determination (POD) using satellite tracking data. 

To achieve the science objectives, especially for measuring the ice-sheet 
topography, the position of the GLAS instrument should be known with an accuracy 
of 5 and 20 cm in radial and horizontal components, respectively. This knowledge 
will be acquired from data collected by the on-board GPS receiver and ground GPS 
receivers and from the ground-based satellite laser ranging (SLR) data. GPS data will 
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be the primary tracking data for the ICESat/GLAS POD, and SLR data will be used 
for POD validation. 


1.2 The POD Problem 

The problem of determining an accurate ephemeris for an orbiting satellite 
involves estimating the position and velocity of the satellite from a sequence of 
observations, which are a function of the satellite position, and velocity. This is 
accomplished by integrating the equations of motion for the satellite from a reference 
epoch to each observation time to produce predicted observations. The predicted 
observations are differenced from the true observations to produce observation 
residuals. The components of the satellite state (satellite position and velocity and 
the estimated force and measurement model parameters) at the reference epoch are 
then adjusted to minimize the observation residuals in a least square sense. Thus, to 
solve the orbit detennination problem, one needs the equations of motion describing 
the forces acting on the satellite, the observation-state relationship describing the 
relation of the observed parameters to the satellite state, and the least squares 
estimation algorithm used to obtain the estimate. 


1.3 GPS-based POD 

Since the earliest concepts, which led to the development of the Global 
Positioning System (GPS), it has been recognized that this system could be used for 
tracking low Earth orbiting satellites. Compared to the conventional ground-based 
tracking systems, such as the satellite laser ranging or Doppler systems, the GPS 
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tracking system has the advantage of providing continuous tracking of a low satellite 
with high precision observations of the satellite motion with a minimal number of 
ground stations. The GPS tracking system for POD consists of a GPS flight receiver, 
a global GPS tracking network, and a ground data processing and control system. 


1.3.1 Historical Perspective 

The GPS tracking system has demonstrated its capability of providing 
high precision POD products through the GPS flight experiment on TOPEX/Poseidon 
(T/P) [Melbourne et al., 1994]. Precise orbits computed from the GPS tracking data 
[ Yu nek et al., 1994; Christensen et al., 1994; Schutz et al., 1994] are estimated to 
have a radial orbit accuracy comparable to or better than the precise orbit 
ephemerides (POE) computed from the combined SLR and DORIS tracking data 
[ Tapley et al., 1994] on T/P. When the reduced-dynamic orbit determination 
technique was employed with the GPS data, which includes process noise 
accelerations that absorb dynamic model errors after fixing all dynamic model 
parameters from the fully dynamic approach, there is evidence to suggest that the 
radial orbit accuracy is better than 3 cm [ Bertiger et al., 1994]. 

While GPS receivers have flown on missions prior to T/P, such as 
Landsat-4 and -5, and Extreme Ultraviolet Explorer, the receivers were single 
frequency and had high level of ionospheric effects relative to the dual frequency T/P 
receiver. In addition, the satellite altitudes were 700 km and 500 km, respectively, 
and the geopotential models available for POD, as they are today, had large errors for 
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such altitudes. As a result, sub-decimeter radial orbit accuracy could not be achieved 
for these satellites. 

Through the GPS flight experiment on T/P several important lessons on 
GPS-based POD have been learned. Those include: 1) GPS Demonstration Receiver 
(GPS/DR) on T/P provides continuous, global, and high precision GPS observable. 
2) GPS-based POD produces T/P radial orbit accuracy similar or better than 
SLR/DORIS. 3) Gravity tuning using GPS measurement was effective [ Tapley et al., 
1996]. 4) Both reduced-dynamic technique and dynamic approach with extensive 
parameterization have been shown to reduce orbit errors caused by mismodeling of 
satellite forces. 


1.3.2 GPS-based POD Strategies 

Several different POD approaches are available using GPS measurements. 
Those include the kinematic or geometric approach, dynamic approach, and the 
reduced-dynamic approach. 

The kinematic or geometric approach does not require the description of 
the dynamics except for possible interpolation between solution points for the user 
satellite, and the orbit solution is referenced to the phase center of the on-board GPS 
antenna instead of the satellite's center of mass. Yunck and Wu [1986] proposed a 
geometric method that uses the continuous record of satellite position changes 
obtained from the GPS carrier phase to smooth the position measurements made with 
pseudorange. This approach assumes the accessibility of P-codes at both the LI and 
L2 frequencies. Byun [1998] developed a kinematic orbit determination algorithm 
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using double- and triple-differenced GPS carrier phase measurements. Kinematic 
solutions are more sensitive to geometrical factors, such as the direction of the GPS 
satellites and the GPS orbit accuracy, and they require the resolution of phase 
ambiguities. 

The dynamic orbit detennination approach [ Tcipley , 1973] requires precise 
models of the forces acting on user satellite. This technique has been applied to many 
successful satellite missions and has become the mainstream POD approach. 
Dynamic model errors are the limiting factor for this technique, such as the 
geopotential model errors and atmospheric drag model errors, depending on the 
dynamic environment of the user satellite. With the continuous, global, and high 
precision GPS tracking data, dynamic model parameters, such as geopotential 
parameters, can be tuned effectively to reduce the effects of dynamic model error in 
the context of dynamic approach. The dense tracking data also allows for the 
frequent estimation of empirical parameters to absorb the effects of unmodeled or 
mismodeled dynamic error. 

The reduced-dynamic approach [Wu et al., 1987] uses both geometric and 
dynamic infonnation and weighs their relative strength by solving for local geometric 
position corrections using a process noise model to absorb dynamic model errors. 

Note that the adopted approach for ICESat/GLAS POD is the dynamic 
approach with gravity tuning and the reduced-dynamic solutions will be used for 
validation of the dynamic solutions. 
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1.4 Outline 

This document describes the algorithms for the precise orbit determination 
(POD) of ICESat/GLAS. Chapter 2 describes the objective for ICESat/GLAS POD 
algorithm. Chapter 3 summarizes the dynamic models, and Chapter 4 describes the 
measurement models for ICESat/GLAS. Chapter 5 describes the least squares 
estimation algorithm and the problem formulation for multi-satellite orbit determination 
problem. Chapter 6 summarizes the implementation considerations for ICESat/GLAS 
POD algorithms. Note that POD ATBD Version 2.2 was written in the pre-launch 
period, and this version (Version 2.3) includes two Appendices to reflect post-launch and 
post-mission POD updates. Note also that contents for Chapter 2 through Chapter 6 are 
the same for Version 2.2 and Version 2.3. Appendix A includes ICESat/GLAS mission 
summary, and the updated POD standards for generating the operational (“Final”) POD. 
Appendix B describes the 2011 POD reprocessing. Bibliography section was updated to 
include references for Appendix A and B. 
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2.0 OBJECTIVE 

The objective of the POD algorithm is to determine an accurate position of 
the center of mass of the spacecraft carrying the GLAS instrument. This position 
must be expressed in an appropriate Earth-fixed reference frame, such as the 
International Earth Rotation Service (IERS) Terrestrial Reference Frame (ITRF), but 
for some applications the position vector must be given in a non-rotating frame, the 
IERS Celestial Reference Frame (ICRF). Thus, the POD algorithm will provide a 
data product that consists of time and the (x, y, z) position (ephemeris) of the 
spacecraft/GLAS center of mass in both the ITRF and the ICRF. The ephemeris will 
be provided at an appropriate time interval, e.g., 30 sec and interpolation algorithms 
will enable detennination of the position at any time to an accuracy comparable to the 
numerical integration accuracy. Furthermore, the transfonnation matrix between 
ICRF and ITRF will be provided from the POD, along with interpolation algorithm. 




3.0 ALGORITHM DESCRIPTION: Orbit 


3.1 ICESat/GLAS Orbit Dynamics Overview 

Mathematical models employed in the equations of motion to describe the 
motion of ICESat/GLAS can be divided into three categories: 1) the gravitational 
forces acting on ICESat/GLAS consist of Earth’s geopotential, solid earth tides, 
ocean tides, planetary third-body perturbations, and relativistic accelerations; 2) the 
non-gravitational forces consist of drag, solar radiation pressure, earth radiation 
pressure, and thermal radiation acceleration; and 3) empirical force models that are 
employed to accommodate unmodeled or mismodeled forces. In this chapter, the 
dynamic models are described along with the time and reference coordinate systems. 


3.2 Equations of Motion, Time and Coordinate Systems 

The equations of motion of a near-Earth satellite can be described in an 
inertial reference frame as follows: 

r =a +a +a (3.2.1) 

where r is the position vector of the center of mass of the satellite, a g is the sum of 
the gravitational forces acting on the satellite, a ng is the sum of the non-gravitational 
forces acting on the surfaces of the satellite, and a is the unmodeled forces which 

act on the satellite due to either a functionally incorrect or incomplete description of 
the various forces acting on the spacecraft or inaccurate values for the constant 
parameters which appear in the force model. 
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3.2.1 Time System 

Several time systems are required for the orbit determination problem. 
From the measurement systems, satellite laser ranging measurements are usually 
time-tagged in UTC (Coordinated Universal Time) and GPS measurements are time- 
tagged in GPS System Time (referred to here as GPS-ST). Although both UTC and 
GPS-ST are based on atomic time standards, UTC is loosely tied to the rotation of the 
Earth through the application of "leap seconds" to keep UT1 and UTC within a 
second. GPS-ST is continuous to avoid complications associated with a 
discontinuous time scale [ Milliken and Zoller, 1978]. Leap seconds are introduced on 
January 1 or July 1, as required. The relation between GPS-ST and UTC is 

GPS-ST = UTC + n (3.2.2) 

where n is the number of leap seconds since January 6, 1980. For example, the 
relation between UTC and GPS-ST in mid-July, 1999, was GPS-ST = UTC + 13 sec. 
The independent variable of the near-Earth satellite equations of motion (Eq. 3.2.1) is 
typically TDT (Terrestrial Dynamical Time), which is an abstract, uniform time scale 
implicitly defined by equations of motion. This time scale is related to the TAI 
(International Atomic Time) by the relation 

TDT = TAI + 32.184 s . (3.2.3) 

The planetary ephemerides are usually given in TDB (Barycentric Dynamical Time) 
scale, which is also an abstract, uniform time scale used as the independent variable 
for the ephemerides of the Moon, Sun, and planets. The transfonnation from the 
TDB time to the TDT time with sufficient accuracy for most application has been 
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given by Moyer [1981]. For a near-Earth application like ICESat/GLAS, it is 
unnecessary to distinguish between TDT and TDB. New time systems are under 
discussion by the International Astronomical Union. This document will be updated 
with these time systems, as appropriate. 


3.2.2 Coordinate System 

The inertial reference system adopted for Eq. 3.2.1 for the dynamic model 
is the ICRF geocentric inertial coordinate system, which is defined by the mean 
equator and vernal equinox at Julian epoch 2000.0. The Jet Propulsion Laboratory 
(JPL) DE-405 planetary ephemeris [Standish, 1998], which is based on the ICRF 
inertial coordinate system, has been adopted for the positions and velocities of the 
planets with the coordinate transformation from barycentric inertial to geocentric 
inertial. 

Tracking station coordinates, atmospheric drag perturbations, and 
gravitational perturbations are usually expressed in the Earth fixed, geocentric, 
rotating system, which can be transformed into the ICRF reference frame by 
considering the precession and nutation of the Earth, its polar motion, and UT1 
transfonnation. The 1976 International Astronomical Union (IAU) precession 
[Lieske et al., 1977; Lieske, 1979] and the 1980 IAU nutation fonnula [ Wahr , 1981b; 
Seidelmann, 1982] with the correction derived from VLBI analysis [Herring et al., 
1991] will be used as the model of precession and nutation of the Earth. Polar motion 
and UT1-TAI variations were derived from Lageos (Laser Geodynamics Satellite) 
laser ranging analysis [ Tapley et al, 1985; Schutz et al, 1988]. Tectonic plate 
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motion for the continental mass on which tracking stations are affixed has been 
modeled based on the AMO-2 model [Minster and Jordan, 1978; DeMets et al., 1990; 
Watkins, 1990]. Yuan [1991] provides additional detailed discussion of time and 
coordinate systems in the satellite orbit determination problem. 


3.3 Gravitational Forces 

The gravitational forces can be expressed as: 

Og = Pgeo + Pst + Pot + Prd + Pn + P,el (3.3.1) 


where 


P = 

r geo 

perturbations due to the geopotential of the Earth 

Pst = 

perturbations due to the solid Earth tides 

Pot = 

perturbations due to the ocean tides 

Pro 

perturbations due to the rotational deformation 

Pn = 

perturbations due to the Sun, Moon and planets 

Prel ~ 

perturbations due to the general relativity 


3.3.1 Geopotential 

The perturbing forces of the satellite due to the gravitational attraction of 

the Earth can be expressed as the gradient of the potential, U, which satisfies the 

2 

Laplace equation, V"U= 0: 
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V U = V (U s + AU st + AU ot + AU r d) = P geo + Pst + Pot + Pra (3.3.2) 
where U s is the potential due to the solid-body mass distribution, A U st is the potential 

change due to solid-body tides, AU ot is the potential change due to the ocean tides, 
and AU rd is the potential change due to the rotational deformations. 

The perturbing potential function for the solid-body mass distribution of 
the Earth, U s , is generally expressed in terms of a spherical harmonic expansion, 
referred to as the geopotential, in a body-fixed reference frame as [Kaula, 1966; 
Heiskanen and Moritz, 1967]: 


oo / 


a. 


TT ( A 3 A GM V- 

U s (r, </>, A) = ^ 2^ 2^ 

r r l=Xm=0 yr J 


PJsm</>)[C lm cos mA + S lm sin m 


iA J 


where 


(3.3.3) 


GM e 

a e 

C/m, S l m 
Plm(smtp) 

r, <f>, A 


the gravitational constant of the Earth 
the mean equatorial radius of the Earth 

nonnalized spherical harmonic coefficients of degree / and order m 
the normalized associated Legendre function of degree / and order 
m 

radial distance from the center of mass of the Earth, the geocentric 
latitude, and the longitude of the satellite 


To ensure that the origin of spherical coordinates coincides with the center of mass of 
the Earth, we define Cio = C\\ =5n =0. 
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3.3.2 Solid Earth Tides 

Since the Earth is a non-rigid elastic body, its mass distribution and the 
shape will be changed under the gravitational attraction of the perturbing bodies, 
especially the Sun and the Moon. The temporal variation of the free space 
geopotential induced from solid Earth tides can be expressed as a change in the 
external geopotential by the following expression [Wahr, 1981a; Dow, 1988; Casotto, 
1989]. 


DM (3) ; 

'M 

l+l 


1+3 

& e 1=2 m=0 k(l,m) 

\ r J 


K r J 



(3.3.4) 


where 


J7W) = (-!)" 


(2/ + 1) — 


]j An ( I + m)\ 


P /m (sin^e 


imX 


f] m (sin (/)) = the unnormalized associated Legendre function of degree / and 


order m 

Hk = the frequency dependent tidal amplitude in meters (provided in 

Cartwright and Tayler [1971] and Cartwright and Edden 
[1973]) 

(9/ ( , Xk = Doodson argument and phase correction for constituent k 

71 

(Xk = 0, if l-m is even; Xk = , if l-m is odd) 

k%, k k = Love numbers for tidal constituent k 

r, (f>, 2 = geocentric body-fixed coordinates of the satellite 
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The summation over k(l,m) means that each different /, m combination has a unique 
list of tidal frequencies, k, to sum over. 

The tidally induced variations in the Earth’s external potential can be 
expressed as variations in the spherical harmonic geopotential coefficients [Eanes et 
al. 1983]. 


AC, 


lm 


A S, 


lm 


(- 1 )"' ^ Tl.Ou J COS 0 A ’ l ~ m even 

a e ^47r(2-S 0m ) T " k i s in ©*• > l~m odd 

(~Y) m x-, n f-sin©,., l-m even 

, 1 ’ Tk° k H k \ k 

a e ^4x(2-S 0m ) k { cosG k , l-m odd 


(3.3.5) 


where 5o m is the Kronecker delta; AQ m and AS are the time -varying geopotential 
coefficients providing the spatial description of the luni-solar tidal effect. 


3.3.3 Ocean Tides 

The oceanic tidal perturbations due to the attraction of the Sun and the 
Moon can be expressed as variations in the spherical harmonic geopotential 
coefficients. The temporal variation of the free space geopotential induced from the 
ocean tide deformation, AU ot , can be expressed as [Eanes et al., 1983] 


oo / 


A U ot = 4 xGp w a e ZZZI 


i+*y- v*‘ 


k 1 = 0 m = 0 


2 / + 1 


a 

V r J 


x l C klm cos (0k±mX) + S' klm sin( 0k±mA)\Pi m (sin^) 


(3.3.6) 
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where p w is the mean density of sea water, k is the ocean tide constituent index, k t is 
the load Love number of degree /, Cz, and Sr, are the unnormalized prograde and 
retrograde tide coefficients, and 0 k is the Doodson argument for constituent k. 

The above variations in the Earth’s external potential due to the ocean tide 
can be expressed as variations in the spherical harmonic geopotential coefficients as 
follows [Eanes et al. 1983], 



k 


^lm } \ B klm 


where F lm , A khn , and B k i m are defined as 




Ana} p w 
M e ■ 


(/+/«)! [ l+kj 

(I-m) ! (2/+ 1 )(2-So, n ) ' 2/+1 


(3.3.7) 


(3.3.8) 


and 


■Akim 

- B k lm 


(C klm C k im ) 

- klm - S klm ) . 


COS 0k + 


(Skim S klm ) 

- klm ~ C klm ) - 


sin 04 - 


(3.3.9) 


3.3.4 Rotational Deformation 

Since the Earth is elastic and includes a significant fluid component, 
changes in the angular velocity vector will produce a variable centrifugal force, 
which consequently deforms the Earth. This deformation, which is called “rotational 
deformation”, can be expressed as the change of the centrifugal potential, U c 
[ Lambeck , 1980] given by 
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U c = ±-a 2 r2 + AU C (3.3.10) 

where 

AU C = z |-(«f+ry?-2^)/ 5 2o(sin^ 

- ^ (coi a>iCosA + oha^sml) P2\{sm.(j)) 

+ j^-[(a^-coi)cos2A - 2ryir«2sin2/l] P22(sin^) (3.3.11) 

and co j = Om\, ah = Onii, or, = and or = (coi+coi+o^). Q is the mean 

angular velocity of the Earth, nij are small dimensionless quantities which are related 
to the polar motion and the Earth rotation parameters by the following expressions: 


m i = x p 


m = -y P 


m 3 


d (UTl-TAI) 
d ( TAI ) 


(3.3.12) 


The first term of Eq. (3.3.10) is negligible in the variation of the 
geopotential, thereby the variation of the free space geopotential outside of the Earth 
due to the rotational deformation can be written as 


AUrd= (^fki AU c (d e ) (3.3.13) 

The above variations in the Earth’s external potential due to the rotational 
defonnation can be expressed as variations in the spherical hannonic geopotential 
coefficients as follows. 

AC 2° = 6GM [ /H|2+;H 2 2 - 2 ( |+ffl 3) 2 ] kTk 2 ~ (\+2m 3 )n 2 k 2 
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Aa -‘ - idk * 3 ^ mn% 

AS21 = ~^rr m 2 (l+m 3 )n 2 k 2 ~ m 2 Q 2 k 2 (3.3.14) 

3 GM e v ' 3 GM e v 

AC 22 = T ^- (m 2 2 -mS-)n 2 k 2 * 0 
AS22 = 777 tt (m 2 mi)n~k 2 « 0 

b(jM e 


As a consequence of Eqs. (3.3.2), (3.3.3), (3.3.4), (3.3.6), and (3.3.13), the 
resultant gravitational potential for the Earth can be expressed as 


t/(r,M) = ^ + ^l£N 

f /=1 ffl=0 V ^ J 


A ( Sin ^ ) 


x [(C/ m +zlC/ m ) cosm/l + (S/ m +AS/ m ) sinm/l] (3.3.15) 

where both the solid Earth and oceans contribute to the periodic variations ZlC/ m and 
AS/ m . 


3. 3. 5 N-Body Perturbation 

The gravitational perturbations of the Sun, Moon and other planets can be 
modeled with sufficient accuracy using point mass approximations. In the geocentric 
inertial coordinate system, the N-body accelerations can be expressed as: 


P n =Y J GM i 


A 

A. 3 


(3.3.16) 


where 
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G = the universal gravitational constant 

= mass of the z'-th perturbing body 

?i = position vector of the z'-th perturbing body in geocentric inertial 

coordinates 

At = position vector of the z'-th perturbing body with respect to the 

satellite 

The values of r, can be obtained from the Jet Propulsion Laboratory Development 
Ephemeris-405 (JPL DE-405) [Standish, 1998]. 


3.3.6 General Relativity 

The general relativistic perturbations on the near-Earth satellite can be 
modeled as [Huang et a/., 1990; Ries et al., 1988], 


Prel 


GM e / 
c-r 3 


(2/3+2 y) 


GM P 


y(? ■ r) 


r + (2+2 y) (? 


r) r ) 


where 


+ 2 (22 x r) 


+ L (\+y) QMe. 


C 2 r 3 LZ' 


(r x r) (?• J) + (7 x J) 


(3.3.17) 


-gm s r es 

2 n3 

C P-ES 

c = the speed of light in the geocentric frame 

7,7 = the geocentric satellite position and velocity vectors 

R es = the position of the Earth with respect to the Sun 


Q 


1 + 2 y 


P-ES X 
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GM e ,GM s = the gravitational constants for the Earth and the Sun, 
respectively 

J = the Earth’s angular momentum per unit mass 

( |/| = 9.8 x 10 8 m 2 /sec) 

L = the Lense-Thirring parameter 

P y = the parameterized post-Newtonian (PPN) parameters 

The first term of Eq. (3.3.17) is the Schwarzschild motion [Huang et ah, 1990] and 
describes the main effect on the satellite orbit with the precession of perigee. The 
second tenn of Eq. (3.3.17) is the effect of geodesic (or de Sitter) precession, which 
results in a precession of the orbit plane [Huang and Ries, 1987]. The last term of 
Eq. (3.3.17) is the Lense-Thirring precession, which is due to the angular momentum 
of the rotating Earth and results in, for example, a 3 1 mas/yr precession in the node of 
the Lageos orbit [Ciufolini, 1986], 


3.4 Nongravitational Forces 

The non- gravitational forces acting on the satellite can be expressed as: 


dig Pdrag^ P solar T P earth P thermal 


where 


Pdrag = perturbations due to the atmospheric drag 

P solar = perturbations due to the solar radiation pressure 

P earth = perturbations due to the Earth radiation pressure 

P thermal = perturbations due to the thermal radiation 


(3.4.1) 
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Since the surface forces depend on the shape and orientation of the satellite, the 
models are satellite dependent. In this section, however, general models are 
described. 


3. 4. 1 Atmospheric Drag 

A near-Earth satellite of arbitrary shape moving with some velocity v in 
an atmosphere of density p will experience both lift and drag forces. The lift forces 
are small compared to the drag forces, which can be modeled as [Schutz and Tapley, 
1980b] 


P(h 


rag 


where 


P 

Vr 

V r 

m 

C d 

A 



= the atmospheric density 
= the satellite velocity relative to the atmosphere 
= the magnitude of iy 
= mass of the satellite 
= the drag coefficient for the satellite 

= the cross-sectional area of the main body perpendicular to vy 


C ,a 

The parameter — — is sometimes referred to as the ballistic coefficient. When more 

1 m 

detailed modeling is needed, the drag force on any specific spacecraft surface, for 
example, the solar panel, can be modeled as 


paneld 


= '2 P 



v r v r 


(3.4.3) 
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where 

Cd p = the drag coefficient for the solar panel 

A p = the solar panel’s area 

y = the angle between the solar panel surface nonnal unit vector, n, 

and satellite velocity vector, v~ r (i.e. cosy= n ■ ) 

\A p cosy\ = the effective solar panel cross sectional area perpendicular to v~ r 
There are a number of empirical atmospheric density models used for 
computing the atmospheric density. These include the Jacchia 71 [. Jacchia , 1971], 
Jacchia 77 [ Jacchia , 1977], the Drag Temperature Model (DTM) [ Barlier et al., 
1977], DTM-2000 [ Bruinsma and Thuillier, 2000], MSIS-90 [ Hedin , 1991] and 
NRLMSISE-00 [ Hedin et al., 1996], The density computed by using any of these 
models could be in error anywhere from 10% to over 200% depending on solar 
activity [ Shum et al., 1986]. To account for the deviations in the computed values of 
density from the true density, the computed values of density, p c , can be modified by 
using empirical parameters which are adjusted in the orbit solution. Once-per- 
revolution density correction parameters [ Elyasberg et al, 1972; Shum et ah, 1986] 
have been shown to be especially effective for these purposes such that 

P = Pc [1 + Ci cos {M+co) + Ci sin (M+<z>)] (3.4.4) 

where 

Ci, Ci = the once-per-revolution density correction coefficients 
M = mean anomaly of the satellite 

co = argument of perigee of the satellite 
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3.4.2 Solar Radiation Pressure 

The Sun emits a nearly constant amount of photons per unit of time. At a 
mean distance of 1 A.U. from the Sun, this radiation pressure is characterized as a 
momentum flux having an average value of 4.56xl0" 6 N lm 2 . The direct solar 
radiation pressure from the Sun on a satellite is modeled as [ Tapley and Ries, 1987] 


P solar P ( 1 T >]) ^ V U 


(3.4.5) 


where 


P = the momentum flux due to the Sun 

77 = reflectivity coefficient of the satellite 

A = the cross-sectional area of the satellite normal to the Sun 

m = mass of the satellite 

v = the eclipse factor ( v = 0 if the satellite is in full shadow, v = 1 if 

the satellite is in full Sun, and 0 < v < 1 if the satellite is in 
partial shadow) 

u = the unit vector pointing from the satellite to the Sun 

Similarly, the solar radiation pressure perturbation on an individual satellite surface, 
like the satellite’s solar panel, can be modeled as 


panels 


= -P V 


|^pCosy| />s ^ 

LJL j l — -\U+ r/p n) 


(3.4.6) 


where 

A p = the solar panel area 

n = the surface nonnal unit vector of the solar panel 
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y = the angle between the solar panel surface nonnal unit vector, n, 

and satellite-Sun unit vector, u (i.e. cos y= u ■ n) 

\ApCOsy | = the effective solar panel cross sectional area perpendicular to u 

The reflectivity coefficient, r), represents the averaged effect over the whole satellite 
rather than the actual surface reflectivity. Conical or cylindrical shadow models for 
the Earth and the lunar shadow are used to determine the eclipse factor, v. Since 
there are discontinuities in the solar radiation perturbation across the shadow 
boundary, numerical integration errors occur for satellites, which are in the 
shadowing region. The modified back differences (MBD) method [. Anderle , 1973] 
can be implemented to account for these errors [Lundberg, 1985; Feulner, 1990]. 


3.4.3 Earth Radiation Pressure 

Not only the direct solar radiation pressure, but also the radiation pressure 
imparted by the energy flux of the Earth should be modeled for the precise orbit 
determination of any near-Earth satellite. The Earth radiation pressure model can be 
summarized as follows [Knocke and Ries, 1987; Knocke, 1989]. 

Pearth = C 1 + ~ Z [( TaE , C0S ^ + eM /<)’’] (3-4.7) 

me j=\ 

where 

rj e = satellite reflectivity for the Earth radiation pressure 

A = the projected, attenuated area of a surface element of the Earth 

A c = the cross sectional area of the satellite 


m 


= the mass of the satellite 
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c = the speed of light 

z = 0 if the center of the element j is in darkness 

1 if the center of the element j is in daylight 
a, e = albedo and emissivity of the element j 

E s = the solar momentum flux density at 1 A.U. 

6 S = the solar zenith angle 

Mg = the exitance of the Earth 

r = the unit vector from the center of the element j to the satellite 

N = the total number of segments 

This model is based on McCarthy and Martin [1977]. 

The nominal albedo and emissivity models can be represented as 

a = ao + a\P\o(sm<f>) + a 2 P 2 o(sm^) (3.4.8) 

e = e 0 + ei.Pio(sin$ + e 2 / 3 2 o(sin^) (3.4.9) 

where 

a\ = co + c\ cos«<f-fo) + ci sinftft-to) (3.4.10) 

e\ = ko + k\ cos<y(Wo) + ^2 simu(t-to) (3.4.11) 

and 

P 10 , P 20 = the first and second degree Legendre polynomial 

(f = the latitude of the center of the element on the Earth’s surface 

co = frequency of the periodic terms (period = 365.25 days) 

/-to = time from the epoch of the periodic term 
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This model, based on analyses of Earth radiation budgets by Stephens et al. [1981], 
characterizes both the latitudinal variation in Earth radiation and the seasonally 
dependent latitudinal asymmetry. 


3.4.4 Thermal Radiation Perturbation 

Since the temperatures of the satellite’s surface are not uniform due to the 
internal and external heat fluxes, there exists a force due to a net thennal radiation 
imbalance. This perturbation depends on the shape, the thermal property, the pattern 
of thennal dumping, the orbit characteristics, and the thermal environment of the 
satellite as a whole. This modeling can be quite complex. For example, if a satellite 
has active louvers for heat dissipation, the thermal force can have specular 
characteristics whereas the heat loss to space from a flat plat is nonnally diffusive. 
Even a clean, perfect spherical satellite like Lageos [Ries, 1989] has been found to 
have a range of detectable thermally induced forces. It is observed for GPS satellites 
that there are unexplained forces in the body-fixed +Y or -Y direction, that is along 
solar panel rotation axis, which causes unmodeled accelerations [Fliegel et al., 1992] 
believed to be of thermal origin. This acceleration is referred to as the “Y-bias”. 
Possible causes of the Y-bias are solar panel axis misalignment, solar sensor 
misalignment, and the heat generated in the GPS satellite body, which is radiated 
preferentially from louvers on the +Y side. Since this Y-bias perturbation is not 
predictable, it can be modeled as 


Pybias & ' ^ Y 


(3.4.12) 
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where uy is a unit vector in the 7-direction, and the scale factor, a, is estimated for 
each GPS satellite. Models, which are satellite-specific, are required to properly 
account for these effects depending on the orbit accuracy needed within a given 
application. 

3. 4. 5 GPS Solar Radiation Pressure Models 

At the 20,000-km altitude of GPS satellite, solar radiation is the dominant 
non-gravitational force acting on the spacecraft. Several GPS solar radiation pressure 
models are currently available, and two of those models are summarized in this 
section. 

Rockwell International Corporation, which was the spacecraft contractor 
for the Block I and II GPS satellites, developed GPS satellite solar radiation pressure 
models, known as ROCK4 for Block I, and ROCK42 for Block II [ Fliegel et al., 
1992]. These models treat a spacecraft as a set of flat or cylindrical surfaces. 
Diffusive and specular forces acting on each surface are computed and summed in the 
spacecraft body-fixed coordinate system. The +Z direction is toward the satellite - 
Earth vector. The +7 direction is along one of the solar panel center beams. The 
satellite is maneuvered so that the X-axis will be kept in the plane defined by the 
Earth, the Sun and the satellite. As a result, the solar radiation pressure forces are 
confined in the X-Z plane, since the 7-axis is perpendicular to the Earth, Sun and the 
satellite plane. The ROCK4 model also provides solar radiation formulas for the X- 
and Z- acceleration components as a function of the angle between the Sun and the 
+Z-axis, e.g. T10 for Block I, and T20 for Block II GPS satellites [. Fliegel et al., 
1992], 
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Recently the Center of Orbit Detennination in Europe (CODE) developed 
a solar radiation pressure (RPR) model by analyzing 5.5 years of GPS orbit solutions 
[Springer et al., 1998], The RPR model is represented by eighteen orbit parameters 
in two different coordinate systems. Those are satellite body-fixed coordinate system 
described above, and the Sun-oriented reference system, which consists of the D-, Y-, 
and /Taxis [ Beutler et al., 1994]. The £)-axis is the satellite-Sun direction positive 
towards the Sun, 7-axis is identical to the ROCK4 7-axis, and /Taxis completes a 
right-handed system. The orbit parameters include three constant terms in the D-, Y-, 
and /^-direction, a once-per-revolution term in the Z-direction, and once- and three 
times-per-revolution terms in the X-direction. The solar radiation acceleration is 
expressed as 

a D = D 0 + D (1 cos (2/3) + D C4 cos(4/3) 
a Y =Y 0 + Y c cos(2(3) 

a B = B 0 + B c cos(2f3) (3.4.13) 

a z = {Z 0 + Z C2 cos(2j3) +Z S2 sin(2/3) 

+Z C4 cos( 4 (3) + Z S4 sin(4 (3)}sin(u - u 0 ) 
a x = {X l0 + X lc cos(2/3) + X [S sin(2(3 )}sin( u - u 0 ) 


+{X 30 + X 3C cos(2/3) + X 3S sin(2 f3)} sin(3u - u 0 ) 
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where u is the argument of latitude of satellite in the orbit plane, no is the latitude of 
the Sun in the orbit plane, and (3 is the angular distance between the orbit plane and 
the Sun. 


3.4. 6 ICESat/GLAS "Box-Wing" Model 

For modeling of non-gravitational perturbations on T/P, the "box-wing" 
model or the so-called macro-model [Marshall et a/., 1992] was developed based on a 
thermal analysis of the spacecraft [Antreasian and Rosborough, 1992], In the macro- 
model, the spacecraft main body and the solar panel are represented by a simple 
geometric model, a box and a wing, and the solar radiation and the thermal forces are 
computed for each surface and summed over the surfaces. For example, the solar 
radiation acceleration for the macro-model is computed using the following equation 
[Milani et al., 1987]. 


P soiar = - p —Tj A i cos0 i 

m 77 


nface 


i = 1 


£ 

2 (y + A cos#, )«, +(! -/?,).? 


(3.4.14) 


where 


P solar 


P 


a 

v 

m 

At 

Oi 


= the solar radiation pressure acceleration 
= the momentum flux due to the Sun 
= the scale factor of the acceleration 
= the eclipse factor (0 for full shadow, 1 for full Sun) 

= mass of the satellite 
= surface area of the i - th plate 

= angle between surface normal and satellite-Sun vector for i - th 
plate 
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n t = surface normal unit vector for /'-th plate 

s = satellite- Sun unit vector 

Pi = specular reflectivity for /'-th plate 

S i = diffusive reflectivity for /'-th plate 

nface = total number of plates in the model 

A similar model is being developed for the ICESat/GLAS satellite, and the model 
parameters, including the specular and diffusive reflectivity coefficients, will be tuned 
using the tracking data. 


3.5 Empirical Forces 

To account for the unmodeled forces, which act on the satellite or for 
incorrect force models, some empirical parameters are customarily introduced in the 
orbit solution. These include the empirical tangential perturbation and the one-cycle- 
per-orbital-revolution (lcpr) force in the radial, transverse, and normal directions 
[Colombo, 1986; Colombo, 1989]. Especially for satellites like ICESat/GLAS which 
are tracked continuously with high precision data, introduction of these parameters 
can significantly reduce orbit errors occurring at the lcpr frequency and in the along 
track direction [Rim etal., 1996]. 


3.5.1 Empirical Tangential Perturbation 

Unmodeled forces in the tangential direction, either along the inertial 
velocity or along the body-fixed velocity, may be estimated by using empirical 
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models during the orbit determination process. This tangential perturbation can be 
modeled empirically as 

Pfangen ~ Ct Mf (3.5.1) 

where 

C t = empirical tangential parameter 

u t = the unit vector in the tangential direction (along inertial velocity 

or body-fixed velocity) 

Such forces are estimated when it is believed that there are mismodeled or unmodeled 
non-conservative forces in the tangential direction. A set of piecewise constants, C t , 
can be estimated to account for these unmodeled tangential perturbations. 


i. 5 . 2 Once-per Revolution R TN Perturbation 

Unmodeled perturbations in the radial, transverse, and normal directions 


can be modeled as 


Pytn 

fa fa 
1 1 

= 


LpJ 



C r cosw + S r sinw 
C t cos u + S t sinw 
C„ cosw + S„ sinw 


where 


(3.5.2) 


Pr 


P, 


Pn 

U 

C r , S r 


= one-cycle-per-revolution radial perturbation 
= one-cycle-per-revolution transverse perturbation 
= one-cycle-per-revolution nonnal perturbation 
= the argument of latitude of the satellite 
= the one-cycle-per-revolution radial parameters 



31 


C t , S t = the one-cycle-per-revolution transverse parameters 
C n , S n = the one-cycle-per-revolution normal parameters 
These empirical perturbations, which are computed in the radial, transverse, and 
nonnal components, are transfonned into the geocentric inertial components. These 
parameters are introduced as needed with complete or subsets of empirical terms 
being used. 
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4.0 ALGORITHM DESCRIPTION: Measurements 

4.1 ICESat/GLAS Measurements Overview 

The GPS measurements will be the primary measurement type for the 
ICESat/GLAS POD, while the laser range measurement will serve as a secondary 
source of verification and evaluation of the GPS-based ICESat/GLAS POD product. 
In this chapter, the mathematical models of the GPS and laser range measurements 
are discussed. 


4.2 GPS Measurement Model 

The GPS measurements are ranges, which are computed from measured 
time or phase differences between received signals and receiver generated signals. 
Since these ranges are biased by satellite and receiver clock errors, they are called 
pseudoranges. In this section, code pseudorange (PR) measurements, phase 
pseudorange measurements (PPR), double-differenced high-low phase pseudorange 
measurements (DDHL) which involve one ground station, two GPS satellites, and 
one low Earth orbiting satellite, are discussed. Consult Hofmann-Wellenhof et al. 
[1992] and Remondi [1984] for more discussion of GPS measurement models. 


4.2. 1 Code Pseudorange Measurement 

The PR measurement, p c pR, can be modeled as follows, 


P pr P ~ c • St t + c- Sty + Spi r0 p + Spi ono + Sp re i 


(4.2.1) 
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where p is the slant range between the GPS satellite and the receiver receiving the 
GPS signal, c is the speed of light, 8t t is the GPS satellite's clock error, 8t r is the 
receiver's clock error, 8p trop is the tropospheric path delay, 8pi ono is the ionospheric 
path delay, and Sp re i is the correction for relativistic effects. 


4.2.2 Phase Pseudorange Measurement 

The carrier phase measurement between a GPS satellite and a ground 
station can be modeled as follows, 

//Ur) = Ah) - m) + mto) (4.2.2a) 

where t Ri is the receive time at the z'-th ground receiver, tj, is the transmit time of the j- 

C j 

th satellite’s phase being received by the z'-th receiver at t Ri , (f> * (t R ) is the computed 
phase difference between the /-th GPS satellite and z'-th ground receiver at t Rp (j) (h) is 
the phase of /-th GPS satellite signal received by z'-th receiver, fy(t R ) is the phase of z'- 
th ground receiver at t Rp to, is the initial epoch of the z'-th receiver, and N/Uo) is the 
integer bias which is unknown and is often referred to as an "ambiguity bias". 
Similarly, the carrier phase measurement between a GPS satellite and a low satellite 
can be modeled as follows, 

<t> C u Ur) = Ah) - UtR) + NJUo) (4.2.2b) 

where t Ru is the received time of the on-board receiver of the user satellite, h u is the 
transmit time of the /-th satellite’s phase being received by the user satellite at t Ru , 
(j) u Ur) is the computed phase difference between j - th GPS satellite and the user 
satellite at t Ru , (f)' U t) is the phase of y'-th GPS satellite signal received by the user 
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satellite, (j) u {tR J is the phase of the user satellite at t Ru , to u is the initial epoch of the 
user satellite, and Nj(to) is the unknown integer bias. 

The signal transmit time of the j- th GPS satellite can be related to the 
signal receive time by 

M = tRi - {piif R yc) - 8tJ. (4.2.3a) 

hi = t Ru - ( pJ(t R ,)/c ) - St'pi (4.2.3b) 

where pj is the geometric line of sight range between j- th GPS satellite and 7-th 
ground receiver, pj is the slant range between j- th GPS satellite and the on-board 
receiver of the user satellite, StJ is the sum of ionospheric delay, tropospheric delay, 

and relativistic effect on the signal traveling from /-th GPS satellite to /-th ground 
receiver, 8t^ is the sum of ionospheric path delay, tropospheric path delay, and 

relativistic effect on the signal traveling from /-th satellite to the on-board receiver of 
the user satellite. Since the time tag, tj or t u , of the measurement is in the receiver 
time scale which has some clock error, the true receive times are 


tR, = ti - 8t Ci (4.2.4a) 

tR u = t u - St Cu (4.2.4b) 

where 8t Ci is the clock error of the /-th ground receiver at t Ri and 8t Cu is the clock error 
of the on-board receiver of the user satellite at t Ru . Since the satellite oscillators and 

the receiver oscillators are highly stable clocks, the (la) change of the frequency over 
the specified period, 4£, is on the order of 10" 12 . With such high stability, the linear 

approximation of <f>(t+ 8t) = <f>(t) + f ■ 8t can be used for 8t which is usually 
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less than 1 second. By substituting Eqs. (4.2.3a) and (4.2.4a) into Eq. (4.2.2a), and 
neglecting higher order terms, Eq. (4.2.2a) becomes 

<fht R ) = An) -/• [ st Ci + P /(t R )/c + ] 

- $(fi) +fi St Ci + Ni(U h ) (4.2.5a) 

Similarly, the phase measurement between the j - th GPS satellite and the user satellite 
can be modeled as follows, 

<l> c J(tR u ) = An) -/• [ 8t Cu + pJ(t Rl )/c + a,! ] 

- UQ +fu 8t Cu + Nj(t 0u ) (4.2.5b) 

By multiplying a negative nominal wavelength, -A = -c //o, where fo is the 

nominal value for both the transmit frequency of the GPS signal and the receiver 
mixing frequency, Eq. (4.2.5a) becomes the phase pseudorange measurement, 

PPR L ! = p/ U r) + j- P/V + j-cdtc-i - 'j^cSt Ci 

-j~- [An) - m] + Cj (4.2.6) 

where Sp^ = c St J and C/ = N/. 

The first tenn of second line of Eq. (4.2.6) can be expanded using the following 
relations: 

An) - m = Ato) - m) + ( if -f,Ut ( 4 . 2 . 7 ) 

J to 

However, p (to) - <Pi(k)) = fi A (to) - fi St Cl (to), which is the time difference 
between the satellite and the receiver clocks at the first data epoch, to. And 
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f 

J to 


[ft - fi) dt is the total number of cycles the two oscillators have drifted apart over 


the interval from to to t,. According to Remondi [1984], this is equivalent to the 
statement that the two clocks have drifted apart, timewise, by amount 

[dtiiti) - St Ci (ti)\ - [stJ(to) - SUtoi Thus, 


Ati) - m = f> • Stc - ■ St Ci (4.2.8) 

After substituting Eq. (4.2.8), Eq. (4.2.6) becomes, 

PPR C / = f p/(t R ) + fjrdpi - £-c Stc + fc 8t Ci + Ci (4.2.9a) 

Jo Jo Jo Jo 

Similarly, the phase pseudorange between y'-th satellite and a user satellite can be 

written as, 

PPR C J fp t i(t Ru ) + f Spj - f c Sti + f c 8t Cu + CJ (4.2.9b) 

To TO TO TO 

Since the GPS satellites have highly stable oscillators, which have 10" 11 or 10" 12 
clock drift rate, the frequencies of those clocks usually stay close to the nominal 
frequency, /q. If the frequencies are expressed as f = fo + A.f J , where Af is clock 
frequency offset from the nominal value, Eqs. (4.2.9a) and (4.2.9b) become as 
follows after ignoring negligible terms: 

PPR C 1 = Pi’UrJ + dpi - c StJ + c St Ci + Cj (4.2.10a) 


PPR C J = pJ(t Ru ) + dpi - c Si + c St Cu + Cj (4.2.10b) 

Note that pj{t R ) and pj(J Ru ) could be expanded as 


p/(t R ) = pj{ti) - pi St,, 


(4.2.11a) 
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PJUrJ ~ PiiUu) ~ Pi! St Cu (4.2.11b) 

Thus, Eqs. (4.2.10a) and (4.2.10b) become 

PPR C / = piiU) + dpi - c dtj + c 8t Ci - p! 8t Ci + C/ (4.2.12a) 

PPR c i! = pm + dpi - c dtc + c dt Cu - p,[ dt Cu + Cj (4.2.12b) 
Eq. (4.2.12a) is the phase pseudorange measurement between a ground receiver and a 
GPS satellite, and Eq. (4.2.12b) is the phase pseudorange measurement between a 
GPS satellite and a user satellite. Note that the clock errors would be estimated for 
each observation epoch. 

4.2. 3 Double-Differenced High-Low Phase Pseudorange Measurement 

By subtracting Eq. (4.2.2b) from Eq. (4.2.2a), a single-differenced high- 
low phase measurement can be formed as follows, 

SDHLPf u = H (t R ) - (4.2.13) 

If another single-differenced high-low phase measurement can be obtained between i- 
th ground receiver, k- th GPS satellite, and the user satellite, a double-differenced 
high-low phase measurement can be formed by subtracting those two single- 
differenced high-low phase measurements. 

DDllLP' t = -f j ■ [dt cl + pjit^lc + dt/_ 

+ /• [ St Cu + pi (t Ru )lc + dti ] 

+ f k ' [ dt Ci + Pi k (t R )/c + dt^f ] 
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- f k - [ 8t Cu + p u k (t Ru )/c + St / ] 

+ /(*,■) - ^*(f f ) - (fi'itu) + </> k (t u ) 

+N* (4.2.14) 

where N J * =N/(t 0i )-Ni(t 0u )-N*(t 0i ) + N k u (t 0u ) . In Eq. (4.2.14), all the phase 

terms associated with the ground station and user satellite receivers are canceled out. 

By multiplying a negative nominal wave length, -A = -c//o, Eq. (4.2.14) 
becomes the double-differenced high-low phase pseudorange measurement, 

ddhcS = Kd • (pKu ) - p‘s \ )) - K-l ■ (W ) - <'«.>) 

V J o ) Wo ) 

Wo ) 

+ c- (W4 

70 / 

+ j£) • < w - W> • (£■) • ( W • W) 

+ Of (4.2.15) 

where Sp^-c- St ^ and C?=-A-N jk . Note that Eq. (4.2.15) contains two 

different time tags, t\ and t u . If the ground station receiver clock and the on-board 
receiver clock are synchronized, then the second line can be canceled out. Since both 
the ICESat/GLAS on-board receiver clock and the ground station receiver clock will 
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be synchronized within 1 microsecond with the GPS System Time, the second line 
can be canceled out. 

Since the GPS satellites have highly stable oscillators, which have 10" 11 or 
10" 12 clock drift rate, the frequencies of those clocks usually stay close to the nominal 
frequency, fo. If the frequencies are expressed as f = fo + Af and f k = fo + Af k , 
Eq. (4.2.15) becomes 

DDHL t = p' i (t li )~ p J u (G ) - p (t Ri ) + p u (t Ru ) 

+ ' Pit (tRu)) - ( Pi k (tR ) - PuUrJ) 

+ £■• ( dl pf-\ <<*«-<*-) 

+ dpi - dpi - dp£ + dp^i 
+ ' ^pi - dpi) - |^— • {dpfr - dpfa ) 

+ Qf (4.2.16) 

For the ICESat/GLAS-GPS case, the single differenced range can be 600 km to 6200 
km. If we assume 10" 11 clock drift rate for GPS satellite clocks, the second line 
contributes an effect, which is at the sub-millimeter level to the double differenced 
range measurement. This effect is less than the noise level, and as a consequence, the 
contribution from the second line can be ignored. Since the perfonnance 
specification of the time-tag errors of the flight and ground receivers for 
ICESat/GLAS mission is required to be less than 1 microsecond with respect to the 
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GPS System Time, the third line also is negligible. The fifth line is totally negligible, 
because even for the propagation delay of 100 m, the contribution from this line is 
less than 10" 9 meters. The first line in Eq. (4.2.16) can be expanded by the linear 
approximation after substituting Eqs. (4.2.4a) and (4.2.4b), to obtain: 

ddhlI = d ft) - ,) - Pi ft ) + P,‘(<, ) 

- [p/ ( ti ) - p/ft )] ■ Stc, + [p/fe) - p/ft)] ■ Stc, 

+ SpJ. - Sp / - op l + .ft/ 

+ C '* (4.2.17) 

This equation is implemented for the double-differenced high-low phase pseudorange 
measurement. The second line does not need to be computed if the ground stations 
and the ICESat/GLAS on-board receiver’s time-tags are corrected in the 
preprocessing stage by using independent clock infonnation from the pseudo-range 
measurement. If such clock information is not available, then the receiver clock 
errors, St Ci and 8t Cu , can be modeled as linear functions, 

Stc, = 8 + bi (U - tio ) (4.2.18a) 

Stc u = a u + b u (t u - t u0 ) (4.2. 1 8b) 

where (a,-, bi) and (a lh b u ) are pairs of clock bias and clock drift for z-th ground station 
receiver clock and the user satellite clock, respectively, and ho and t u o are the 
reference time for clock parameters for z-th ground station receiver clock and the user 
satellite clock. 
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The third line of Eq. (4.2.17) includes the propagation delay and the 
relativistic effects for the high-low phase converted measurement. These effects are 
discussed in more detail in the following sections. 


4. 2. 4 Corrections 

4. 2. 4. 1 Propagation Delay 

When a radio wave is traveling through the atmosphere of the Earth, it 
experiences a delay due to the propagation refraction. Atmospheric scientists usually 
divide the atmosphere into four layers: the troposphere, the stratosphere, the 
mesosphere, and the thennosphere. The troposphere, the lowest layer of the Earth’s 
atmosphere, contains 99% of the atmosphere’s water vapor and 90% of the air mass. 
The tropospheric bending is therefore treated using both dry and wet components. 
The dry path delay is caused by the atmosphere gas content along the propagated path 
through the troposphere while the wet path delay is caused by the water vapor content 
along the same path. Since the tropospheric path delay of a radio wave is frequency 
independent, this path delay cannot be isolated using multiple frequencies. The 
tropospheric path delays caused by the dry portion, which accounts for 80% or more 
of the delay, can be modeled with an accuracy of two to five percent for L-band 
frequencies [ Atlshuler & Kalaghan, 1974], Although the contribution from the wet 
component is relatively small, it is more difficult to model because surface 
measurements of water vapor cannot be applied to completely describe the regional 
variations in the water vapor distribution, especially with respect to horizontal 



42 


variation, of the water vapor field. There are several approaches to model the wet 
component of the tropospheric path delay. One approach is to use one of the 
empirical atmospheric models based on the measurement of meteorological 
parameters at the Earth’s surface or the altitude profile with radisondes and apply 
regional modeling. The other approach is to map the water vapor content in various 
directions directly using devices like water vapor radiometer (WVR). List of 
references for these approaches can be found in Tralli et al. [1988]. A third approach 
is to solve for tropospheric path delay parameters. Chao’s model [Chao, 1974], 
modified Hopfield model [ Goad and Goodman, 1974; Remondi, 1984], or MTT 
model [Herring, 1992] are among several candidates which can be implemented for 
the tropospheric correction. 

The ionosphere is a region of the Earth’s upper atmosphere, approximately 
100 km to 1000 km above the Earth’s surface, where electrons and ions are present in 
quantities sufficient to affect the propagation of radio waves. The path delay will be 
proportional to the number of electrons along the slant path between the satellite and 
the receiver, and the electron density distribution varies with altitude, time of day 
time of year, solar and geomagnetic activity, and the time within the solar sunspot 
cycle. The ionospheric path delay depends on the frequency of the radio signal. The 
ionospheric bending on LI GPS measurement will vary from about 0.15 m to 50 m 
[Clynch and Coco, 1986]. Some of this delay can be eliminated by ionospheric 
modeling [for example, Finn and Matthewman, 1989]. However, more accurate 
corrections can be made by using the dual frequency measurements routinely 
acquired by the GPS receivers. The correction method for the dual frequency GPS 
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measurements can be found in Section 6.5.2. Hofmann-Wellenhof et al. [1992] 
provides more detailed description of the propagation delay for GPS measurements. 


4.2. 4.2 Relativistic Effect 

The relativistic effects on GPS measurements can be summarized as 
follows. Due to the difference in the gravitational potential, the satellite clock tends 
to run faster than the ground station’s [ Spilker , 1978; Gibson, 1983], These effects 
can be divided into two parts: a constant drift and a periodic effect. The constant drift 
can be removed by off-setting the GPS clock frequency low before launch to account 
for that constant drift. The periodic relativistic effects can be modeled for a high-low 
measurement as 


Apsrel =%(rr v i- n, ■ n) (4.2.20) 

where 

Ap S rei = correction for the special relativity 
c = speed of light 

?/, v/ = the position and velocity of the low satellite or tracking stations 

Vh = the position and velocity of the high satellite 

The coordinate speed of light is reduced when light passes near a massive body 
causing a time delay, which can be modeled as [Holdridge, 1967] 


Apgrel 


= (1 + r)^h In 


r tr T r rec "P P ) 
r tr r rec ~ P I 


(4.2.21) 


where 


= correction for the general relativity 
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y = the parameterized post-Newtonian (PPN) parameter (y = I for 

general relativity) 

GM e = gravitational constant for the Earth 

p = the relativistically uncorrected range between the transmitter 

and the receiver 

r tr = the geocentric radial distance of the transmitter 

r rec = the geocentric radial distance of the receiver 


4.2. 4.3 Phase Center Offset 

The geometric offset between the transmitter and receiver phase centers 
and the effective satellite body-fixed reference point can be modeled depending on 
the satellite orientation (attitude) and spacecraft geometry. The ICESat/GLAS 
antenna location will be known and implemented when the fabrication of the satellite 
is complete. However, the location of the antenna phase center with respect to the 
spacecraft center of mass will also be required. This position vector will be 
essentially constant in spacecraft fixed axes, but this correction is necessary since the 
equations of motion refer to the spacecraft center of mass. 


4.2. 4.4 Ground Station Related Effects 

In computing the double-differenced phase-converted high-low pseudo- 
range measurement, it is necessary to consider the effects of the displacement of the 
ground station location caused by the crustral motions. Among these effects, tidal 
effects and tectonic plate motion effects are most prominent. 
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Station displacements arising from tidal effects can be divided into three 

parts, 

Atide ~ ^dtide ^ ocean T A rotate (4.2.22) 

where 

A tide = the total displacement due to the tidal effects 

A dude = the displacement due to the solid Earth tide 

A ocean = the displacement due to the ocean loading 

A m tate ~ the displacement due to the rotational deformation 

The approach of the IERS Conventions [McCarthy, 1996] have been implemented for 
the solid Earth tide correction. Ocean loading effects are due to the elastic response 
of the Earth’s crust to loading induced by the ocean tides. The displacement due to 
the rotational deformation is the displacement of the ground station by the elastic 
response of the Earth’s crust to shifts in the spin axis orientation [Goad, 1980] which 
occur at both tidal and non-tidal periods. Detailed models for the effects of solid 
Earth tide, the ocean loading, and the rotational deformation, can be found in Yuan 
[1991], 

The effect of the tectonic plate motion, which is based on the relative plate 
motion model AMO-2 of Minster and Jordan [1978], is modeled as 

A „=Kx\)((-g (4.2.23) 

where 

A tect = the displacement due to the tectonic motion 
c o p = the angular velocity of the tectonic plate 

R So = the Earth-fixed coordinates of the station at /, 
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to = a reference epoch 


4.2. 5 Measurement Model Partial Derivatives 

The partial derivatives of Eq. (4.2.18) with respect to various model 
parameters are given in this section. The considered parameters include the ground 
station positions, GPS satellite’s positions, ICESat's positions, clock parameters, 
ambiguity parameters, and tropospheric refraction parameters. 


The partial derivatives of Eq. (4.2.18) with respect to the z'-th ground 
station positions, (xu, X 2 i, x 3 ,), are 


dDDHL 

8X mi 


■jk 


(x mi -x m J ) (x mi -x m ) 


Pi 


Pi 


for m=l ,2,3 (4.2.24) 


where pi is the range between z'-th ground station receiver and y'-th transmitter, and 
p k is the range between z'-th ground station receiver and k - th transmitter such that 


Pi = V(x h - - x\i) 2 + (x 21 - xP) 2 + (x 3i - xpf 

(4.2.25) 

Pi k = V (x\i -Xi k ) 2 + (x 2 i - X 2 k ) 2 + (x 3i - x 3 k ) 2 

(4.2.26) 


and (x ]/, xP, xp) and (x \ k , xi k , X 3 k ) are the j - th and k - th transmitter Cartesian 
positions, respectively. 


The partial derivatives of Eq. (4.2.18) with respect to the y-th and k - th 


transmitter positions are 
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J k 


dDDHIftu _(x mi -x m ) , (x mu -x m ) 


cx„ 


Pi 


Pu 


SDDHL‘ t (x„,-x„‘) (x r „,-Q 

k k 

Pi 


k ’ 


for m=l,2,3 (4.2.27) 


for m=l,2,3 (4.2.28) 


ox m pi p;, 

where pj is the range between j - th transmitter and the user satellite, and pj' is the 
range between A:-th transmitter and the user satellite such that 

Pli = V (xiu -Xi jf + ( x lu - Xpf + (x 3u - X 3 J ) 2 (4.2.29) 

(4.2.30) 


Pu = V(x hi - Xi k ) 2 + (x lu - x 2 k f + (x 3u - x 3 k f 
and (xi u , X 2 u, xt, u ) are the user satellite's Cartesian positions. 


The partial derivatives of Eq. (4.2.18) with respect to the user satellite 
positions are 


dDDHL 

dx.„„ 


Jk 




Pu 


Pu 


for m=l,2,3 (4.2.31) 


The partial derivatives of Eq. (4.2.18) with respect to the clock parameters 
of Eqs. (4.2.19a) and (4.2.19b) are 


dDDIIL 

da 


dPPHL 

db~ 


Jk 


Jk 


=~(P/-Pi) 

(4.2.32) 

= -(p/-PiHti-tio) 

(4.2.33) 


and 
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dDDHL 

da„ 


Jk 


= {pi-pt) 


dDDHL 1 

dh~ 


Jk 


■ = (pl-p k «)-(t«- t u o) 


(4.2.34) 


(4.2.35) 


The partial derivative of Eq. (4.2.18) for the double-differenced ambiguity 
parameter, C ]k , is 


dDDHDh, 

dC Jk 

ill 


(4.2.36) 


When Chao’s model is used, the partial derivative of Eq. (4.2.18) with 
respect to the z'-th ground station’s zenith delay parameter, Z,, is 

f \ 


dDDHL 

azT 


Jk 


1 




1 


. , 0.00143 . 0.00035 

sinE, -l : sinE, H : 

tanii/ +0.0445 tanE/ +0.017 


. 0.00143 

sinC. + 


■ + 


tan E +0.0445 


. k 0.00035 
sinC + 


tan E. +0.017 


(4.2.37) 

where Ej and E-/ ( are the elevation angles of the y-th and k - th GPS satellite 
transmitters from z'-th ground station, respectively. 
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4.3 SLR Measurement Model 

4.3.1 Range Model and Corrections 

Laser tracking instruments record the travel time of a short laser pulse 
from the reference point (opticalaxis) to the satellite retroreflector and back. The 
one-way range from the reference point of the ranging instrument to the retroreflector 
of the satellite, p°, can be expressed in terms of the round trip light time, A i as 

p° = -cAt + £ (4.3.1) 

where 

c = the speed of light 

a = measurement error. 

The computed one-way signal path between the reference point on the 
satellite and the ground station, p c , can be expressed as 

p — |r - r s .| + Ap lrop + Ap^ + Ap c m (4.3.2) 

where 

r = the satellite position in geocentric coordinates 

r s = the position of the tracking station in geocentric coordinates 

Aptrop = correction for tropospheric delay 

Apgrei = correction for the general relativity 

Apc.m. = correction for the offset of the satellite's center-of-mass and the 

laser retroreflector 

The tropospheric refraction correction is computed using the model of Marini and 
Murray [1973]. The correction for the general relativity in SLR measurements is the 
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same as for GPS measurement, which is expressed in Eq. (4.2.21). The effects of the 
displacement of the ground station location caused by the crustral motions should be 
considered. These crustral motions include tidal effects and tectonic plate motion 
effects, which are described in Eqs. (4.2.22) and (4.2.23), respectively. 


4. 3. 2 Measurement Model Partial Derivatives 

The partial derivatives of Eq. (4.3.2) with respect to various model 
parameters are derived in this section. The considered parameters include the ground 
station positions, satellite’s positions. 


The partial derivatives of Eq. (4.3.2) with respect to the ground station 
positions, (r sU r s2 , r s3 ), are 


d P C _ (G ~ r i) 
Sp P 


for i=l, 2, 3 (4.3.3) 


where (r \ , r 2 , r 3 ) are the satellite's positions, and p is the range between the ground 
station and the satellite such that 


P = V (n - r s i) 2 + (r 2 - r s2 ) 2 + (r 3 - r s3 ) 2 


(4.3.4) 


The partial derivatives of Eq. (4.3.2) with respect to the satellite's 
positions, (r \ , r 2 , r 3 ), are 


dp c _ Q; -r„.) 
dr. 


P 


for /= 1 ,2,3 


(4.3.5) 
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5.0 ALGORITHM DESCRIPTION: Estimation 

A least squares batch filter [Tap ley, 1973] is our adopted approach for the 
estimation procedure. Since multi-satellite orbit determination problems require 
extensive usage of computer memory for computation, it is essential to consider the 
computational efficiency in the problem formulation. This section describes the 
estimation procedures for ICESat/GLAS POD, including the problem formulation for 
multi-satellite orbit determination. 

5.1 Least Squares Estimation 

The equations of motion for the satellite can be expressed as 

Xit) = F(X,t), X(to) = X 0 (5.1.1) 

where X is the //-dimensional state vector, A is a non-linear /7-dimensional vector 
function of the state, and Xq is the value of the state at the initial time to, which is not 
known perfectly. The tracking observations can be expressed as discrete 
measurements of quantities, which are a function of the state. Thus the observation- 
state relationship can be written as 

Y i =G(X i ,t i )+ % i = (5.1.2) 

where Y t is a p vector of the observations made at time /„ (A/, tj) is a non-linear vector 
function relating the state to the observations, and Si is the measurement noise. 

If a reference trajectory is available and if X, the true trajectory, and X* , 
the reference trajectory, remain sufficiently close throughout the time interval of 
interest, the trajectory for the actual motion can be expanded in a Taylor series about 



52 


the reference trajectory to obtain a set of differential equations with time dependent 
coefficients. Using a similar procedure to expand the nonlinear observation-state 
relation, a linear relation between the observation deviation and the state deviation 
can be obtained. Then, the nonlinear orbit determination problem can be replaced by 
a linear orbit determination problem in which the deviation from the reference 
trajectory is to be determined. In practice, this linearization of the problem requires 
an iterative adjustment which yields successively smaller adjustments to the state 
parameters to optimally fit the observations. 

Let 


x(t) = X(t)-X*(t) y(t)= Y(t)-Y*(t) (5.1.3) 

where X* (t) is a specified reference trajectory and Y*(t) is the value of the observation 
calculated by using X* ( t ). Then, substituting Eq. (5.1.3) into Eqs. (5.1.1) and (5.1.2), 
expanding in a Taylor's series, and neglecting higher order terms leads to the relations 


x = A(t)x, 


x(to ) = x 0 


yi = Hpa+ St 




where 


A(t) = ^{X*, f) 
dX 


~ dG * 

H ( X , t) 

dX 


The general solution to Eq. (5.1.4) can be expressed as 


x(t) = <I(t,t 0 )xo 


(5.1.4) 


(5.1.5) 


(5.1.6) 


where the state transition matrix cS(t,to) satisfies the differential equation: 
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WJo) = A(t)<Ht,to), Wo, t 0 ) = 1 (5.1.7) 

where / is the n xn identity matrix. 

Using Eq. (5.1.5), the second of Eq. (5.1.3) may be written in terms of the 
state at to as 

Yi = HiWiJotyo + (5.1.8) 

Using the solution for the linearized state equation (Eq. (5.1.6)), Eq. (5.1.8) may be 
rewritten as 

y=Hxo + a (5.1.9) 

where 

"ul H^{t v t 0 ) U 

y = : H= i i (5.1.10) 

_y i _ H^{ti,t 0 ) s t 

where y and a are in vectors (m = Ixp ) and H is an mxn matrix. Equation (5.1.9) is a 
system of m equations in n unknowns. In practical orbit determination problems, 
there are more observations than estimated parameters (m > n), which means that Eq. 
(5.1.9) is overdetennined. It is usually assumed that the observation error vector, a, 
satisfies the a priori statistics, E[a\ = 0 and E[ss T ~\ = W~ l . By scaling each term in Eq. 
(5.1.9) by W 1/2 , the condition 

j^l/2 ^T-jfpT/2 = yy\!2 yy-\ yyT! 2 = j 


is obtained. 


(5.1.11) 
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An approach to obtain the best estimate of x, given the linear observation- 
state relations (Eq. (5.1.9)) is described in the following discussions. The method 
obtains the solution by applying successive orthogonal transfonnations to the linear 
equations given in Eq. (5.1.9). Consider the quadratic performance index 

J=\ I W m {Hx -y)\ 2 =±(Hx -y) T W(Hx - y) (5.1.12) 

The solution to the weighted least-squares estimation problem (which is 
equivalent to the minimum variance and the maximum likelihood estimation problem, 
under certain restrictions) is obtained by finding the value x which minimizes Eq. 
(5.1.12). To achieve the minimum value of Eq. (5.1.12) let Q be an mxm orthogonal 


matrix. Hence, it follows that Eq. (5.1.12) can be expressed as 

J=\\QW m {Hx- y) \- (5.1.13) 

Now, if Q is selected such that 

QW m H= ^ QW m y= b (5.1.14) 

where R is nxn upper-triangular, 0 is an (in-n)xn null matrix, b is nx 1 vector, and e is 
an ( m-n)x 1 vector. Equation (5.1.13) can be written then as 

J(x) = ±\Rx-b\ 2 +±\e[ 2 (5.1.15) 

The value ofx, which minimizes Eq. (5.1.12), is obtained by the solution 

Rx = b (5.1.16) 

and the minimum value of the performance index becomes 

J(x) = L\ e \2=L\y_HSY (5.1.17) 

That is, e provides an estimate of the residual error vector. 
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The procedures are direct and for implementation requires only that a 
convenient computational procedure for computing QW l,2 H and QW l,2 y be 
available. The two most frequently applied methods are the Givens method, based on 
a sequence of orthogonal rotations, and the Householder method, based on a series of 
orthogonal reflections [Lawson and Hanson, 1974]. 

In addition to the expression for computing the estimate, the statistical 
properties of the error in the estimate, R , are required. If the error in the estimate, /], 
is defined as 


rj = x-x (5.1.18) 

it follows that 

E[rj\ = E[x-x\ = E[R' l b-x] (5.1.19) 

Since 


QW m y = QW m Hx + QW m e 


leads to 


b = Rx + e (5.1.20) 

it follows that 

E[rj~\ = E[R l (Rx + s)- x] = E[R l s~\ (5.1.21) 

As noted in Eq. (5.1.11), if the observation error, a, is unbiased, s= QW l/1 a will be 
unbiased and 


E[rj\ = 0 


(5.1.22) 
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Hence, x will be an unbiased estimate of x. Similarly, the covariance matrix for the 
error in x can be expressed as 

P = E[t]r/ r ] 

= E[R - 1 ss T R - T ] = R - l E[ss T ]R - t (5.1 .23) 

If the observation error, e, has a statistical covariance defined as E\ss T ~\ = W' 1 , the 
estimation error covariance matrix is given by 

E[ss T ] = W l/2 E[s£ T ]W T/2 = IV 12 W~ l W 72 = 1. Consequently, relation (5.1.23) leads 
to 


P = R a R- t (5.1.24) 

It follows then that the estimate of the state and the associated error covariance matrix 
are given by the expressions 

x = R A b (5.1.25) 

P = R a R- t (5.1.26) 


5.2 Problem Formulation for Multi-Satellite Orbit Determination 

Proper categorization of the parameters will help to clarify the problem 
formulation. Parameters can be divided into two groups: dynamic parameters and 
kinematic parameters. Dynamic parameters need to be mapped into other states by 
using the state transition matrix, which is usually computed by numerical integration, 
while kinematic parameters are treated as constant throughout the computation. 
Dynamic parameters can be grouped again into two parts as the local dynamic 
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parameters and global dynamic parameters. Local dynamic parameters are satellite- 
specific. Global dynamic parameters are parameters, which influence every satellite, 
such as those defining gravitational forces. 

Following the categorization described above, the estimation state vector 
is defined as 


Xkp 
X ss 
Xldp 
Xgdp 


where 


(5.2.1) 


Xkp = the kinematic parameters (npp ) 

Xss = the satellite states (n ss ) 

Xldp = the local dynamic parameters (np/p) 

Xqdl = the global dynamic parameters (n gc / p ) 
and X ss consists of satellites’ positions and velocities, i.e. X ss = \Xpos-> XvelY • For 


ns-satellites, where ns is the total number of satellites which will be estimated, Xss 



becomes 
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where r, and v; are the 3x1 position and velocity vectors of the i-th satellite, 
respectively. 

The differential equations of state, Eq. (5.1.1), becomes 


m = F(X,t) 


' 0 ' 

Xss 

0 ’ 

L 0 J 


X(to) = X 0 


(5.2.2) 


where 


X.. = 


f 



(5.2.3) 


and fi = a gi + a ngl for /'-th satellite. Eq. (5.2.2) represent a system of n nonlinear first 
order differential equations which includes n ss = 6 xns of Eq. (5.2.3). After the 
linearization process described in section 5.1, Eq. (5.2.2) becomes Eqs. (5.1.6) and 
(5.1.7). 

Since Eq. (5.1.7) represents n coupled first order ordinary differential equations, the 
dimension of the integration vector becomes n ss + n . However, A(t) matrix is a 
sparse matrix, because of the nature of the parameters. And A(t) matrix becomes 
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even sparser, since each satellite’s state is independent of the others, i.e. ( 7 /, vf) is 

independent of (?/, vj) for i *j. Using the partitioning of Eq. ( 5 . 2 . 1 ), Aft) becomes 

'0 0 0 0 O' 

0 0/00 

A = 0 A 32 A 33 A 34 A 35 ( 5 . 2 . 4 ) 

o 0 0 0 0 

L0 0 0 0 0 J 

where 



Note that A 32, A 33, and A 34 are all block diagonal matrix, and A 33 would be zero if 
the perturbations do not depend on satellites’ velocity. Atmospheric drag is one 
example of perturbations, which depend on the satellite’s velocity. 

If & = for i,j = 1 , • 5 , Eq. ( 5 . 1 . 7 ) becomes 
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' 0 0 0 0 O' 

031 032 033 034 035 

0 = Bn B\2 B 13 B l4 B 15 (5.2.5) 

0 0 0 0 0 

L o o o o o J 

where By = A 32 <hj + ^3303/ + -4 34 0 4y - + ^ 35 05y fory = 1, , 5. 

Integrating the first row and last two rows of Eq. (5.2.4) with the initial 
conditions, (IK to, to) = I yields the results that 0n=0 44 = 055=4 and 

01 2 = 01 3 = 01 4 = 01 5 = 04 1 = 042 = 043 = 045 = 05 1 = 052 = 053 = 054 =0 . After substituting these 
results to By,j =1, • -,5, we have 

4?11 = -432021 + -4 33 05i 

B 12 = -4 32 022 + -4 33 032 

Bn = -4 32 023 + -4 33 053 (5.2.6) 

#14 = -4 32 024 + -4 33 03 4 + -4 34 
4? 15 = -4 32 025 + -4 33 035 + -4 35 

From Eq. (5.2.5) and Eq. (5.2.6), we have 

021 = 031 (5.2.7a) 


022 - 032 


(5.2.8a) 



61 


023 = hi (5.2.9a) 

02 4 = 034 (5.2.10a) 

025 = hi (5.2.11a) 

<h\ = -432021 + -433 031 (5.2.7b) 

hi = -4 32 hi + -4 33 hi (5.2.8b) 

hi = -4 32 hi + Aahi (5.2.9b) 

034 = -4 32 024 + -433 034 + ^34 (5.2.10b) 

035 = -432025 + Aiihi + -435 (5.2.11b) 

From Eqs. (5.2.7a) and (5.2.7b) 

021 --433 021 --432 021 = 0, 02l(O) = O 021 (0) = 0 (5.2.12) 

If we define the partials of accelerations with respect to each group of parameters for 
the /-th satellite as follows, 

— = (5.2.13a) 

d?i 

— = DAD Vi (5.2.13b) 

3 ?/ 
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Vi 

dX LDP, 


= DLDPj 


dfi 

dX GDP, 


= DGDPj 


(5.2.13c) 


(5.2.13d) 


and <jh\ is partitioned as (jh\ =[^i 15 • • •, 0>i by 3 x n kp submatrix, (fh\ p then, Eq. 
(5.2.12) become 

021,- - DADVi fa. - DADRj fay. = 0, i = 1, ■ -,ns (5.2.14) 


After applying the initial conditions, 02i,(O) = 0 and 02i,(O) = 0, to Eq. (5.2.14), we 
have 02i = 0. And from Eq. (5.2.7a) = 0. From Eqs. (5.2.8a) and (5.2.8b), and 

Eqs. (5.2.9a) and (5.2.9b), we have similar results as follows. 

022, - DADVi 022, - DADRj ;2 . = 0, i= 1, ■ -,ns (5.2.15) 

023, - DADVj 023, - DADRj 02 3i =0, i= 1 , • -,ns (5.2.16) 

with the initial conditions 022,(0) = 1 , 022,(0) = 0, 023 ,(0) = 0, and 023 ,(0) = / for i = 
1 , ■ -,ns. 

From Eqs. (5.2.10a) and (5.2.10b), we have 
024 ' ^33 024 - ^32 024 = ^34, 


024(0) = 0 024(0) = 0 (5.2.17) 
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r ~\ T 

If 024 is partitioned as 0 24 = 0 24 , ■ ■ • , 0 24 with 3 x n ldp submatrix, where n ldp is 

the z'-th satellite’s number of local dynamic parameters, then it can be shown that all 
the off-block diagonal terms become zero and the above equation becomes, 

024,. - DAD Vi 024, - DADRi 024, = DLDP t , i= 1 ,---,ns (5.2.18) 

with the initial conditions 024,(0) = 0 and 024,(0) = 0 for i = 1, • • • ,ns. 

From Eqs. (5.2. 1 la) and (5.2. 1 lb), we have similar results for 025. 

025, - DADVi 025, - DADRi 025, = DGDP U /'= 1, ■ • • ,ns (5.2.19) 

with the initial conditions 025,(0) = 0 and 025,(0) = 0 for i = 1, • • • ,ns. 

Combining all these results, we have the state transition matrix for multi- 
satellite problem as follows: 

^ 0 0 0 0 0 ^ 

021 022 023 024 025 

® — 02j 0,2 023 024 025 (5.2.20) 

0 0 0 0 0 

v 0 0 0 0 o y 

where 0 21 = 0 21 = 0 and 
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024 1 

0 


025 1 



025 - 


0 

024, K . 


. 025™. 


By defining 0: and 07, for 7-th satellite as follows, 

0i = [022, 023, 024 i 025] (5.2.21a) 

0v, = [022 , 023 , 024 , 025 ] (5.2.21b) 


we can compute 07,. = [ 022 , 023, 024, 02 sJ by substituting Eqs. (5.2. 15)-(5.2. 16) 
and Eqs. (5.2.1 8)-(5.2. 19). 


DAD Vi 022 i + DADRj 022, 


0E 


DADVi 023, + DA DR, 0 23 ,. 
DAD Vi 024,. + DADRj 024,.+ DLDPj 


DAD Vi 0,5 , + DA DR, 0> 5 ,. + DGDP, 


(5.2.22) 


After rearranging this equation, we get 

07, = DADVi 07, + DADRj 0, + [0 3 . y3 0 3x3 DLDPj DGDP,] (5.2.23) 


Eq. (5.2.23) represents 3 x (6+ni ( / l ,+n g( / p ) first order differential equations for the 7-th 
satellite. Therefore, the total number of equations for ns satellites 

ns 

becomes X 3x (6 + n ldp) + n gdpi ) . 

i= 1 
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Since multi-satellite orbit determination problem includes different types 
of satellites in terms of their perturbations and integration step size, a class of satellite 
is defined as a group of satellites which will use the same size of geopotential 
perturbation and the same integration order and step size. For /-classes of satellites, 
the integration vector, X is defined as 
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XlNT 


n i 

n ns i 

Ml 

ns i 
VJ 1 

V 1 ns i 

Ml 
M ns i 

r/1 

r lnsi 

M 

Ms, 

V/l 

v Insi 

Mi 

^V/ns, 


(5.2.24) 


where ns,- is the number of satellites for 7-th class, ry and v,y are the position and 
velocity of the y'-th satellite of 7-th class, respectively. <fc is the state transition 
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matrix for the y'-th satellite’s positions of 7-th class and <fc.. is the state transition 


matrix for the y-th satellite’s velocities of 7-th class. 



(5.2.25) 
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Eq. (5.2.25) is numerically integrated using a procedure such as the Krogh-Shampine- 
Gordon fixed-step fixed-order formulation for second-order differential equations 
[Lundberg, 1981] for each class of satellites. For the ICESat/GLAS-GPS case, two 
classes of satellites need to be defined. One is for the high satellites, e.g. GPS, and 
the other is for the low satellite, e.g. ICESat/GLAS. 


5.3 Output 

Although a large number of parameters are available from the estimation 
process as given by Eq. (5.2.24), the primary data product required for the generation 
of other products is the ephemeris of the ICESat/GLAS spacecraft center of mass. 
This ephemeris will be generated at a specified interval, e.g., 30-sec and will include 
the following: 

t in GPS time 

3 position components of the spacecraft center of mass in ICRF and ITRF 
7)™/ the 3x3 transformation matrix between ICRF and the ITRF. 

The output quantities will be required at times other than those contained in the 
generated ephemeris file. Interpolation methods, such as those examined by 
Engelkemier [1992] provide the accuracy comparable to the numerical integration 
accuracy itself. With these parameters the ITRF position vector can be obtained as 
well by forming the product of the transformation matrix and the position vector in 


ICRF. 
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6.0 IMPLEMENTATION CONSIDERATIONS 

In this chapter, some considerations for implementing ICESat/GLAS POD 
algorithms are discussed. Section 6.1 describes the POD software system in which 
the POD algorithms are implemented, and the necessary input files for the software 
are defined. Section 6.2 describes the POD products. Section 6.3 describes the 
ICESat/GLAS orbit and attitude. Section 6.4 discusses the expected ICESat/GLAS 
orbit accuracy based on simulations. Section 6.5 summarizes the POD processing 
strategies. Section 6.6 discusses the plans for pre-launch and post-launch POD 
activities. Section 6.7 considers computational aspects. 


6.1 POD Software System 

The POD algorithms described in the previous chapters were implemented 
in a software system, referred to as MSODP1 (Multi-Satellite Orbit Determination 
Program 1). This software has been developed by the Center for Space Research 
(CSR), and shares heritage with UTOPIA [ Schutz and Tapley, 1980a]. This software 
can process SLR data and Doppler data in addition to GPS pseudo-range and double- 
differenced carrier phase data. A version of this POD software will be placed under 
change control at ICESat/GLAS launch. MSODP1 requires input files, some of 
which define model parameters, and the following section discusses these necessary 
input files. 
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6. 1. 1 Ancillary Inputs 

Some model parameters require continual updating through acquisition of 
input information hosted on various standard anonymous ftp sites. This includes the 
Earth orientation parameters, x p , y p , and UT1, and solar flux data. Other files, which 
are considered to be static once "tuned" to ICESat/GLAS requirements include the 
planetary ephemerides, geopotential parameters, and ocean tides parameters,. In 
addition, information about the spacecraft attitude is required for the box-wing 
spacecraft model in the computation of non-gravitational forces and to provide the 
correction for the GPS phase center location with respect to the spacecraft center of 
mass. The real-time attitude obtained during flight operations is thought to be 
adequate for this purpose, but it will be checked against the precise attitude during the 
Verification Phase. Also, the GPS data from the IGS ground network and the 
ICESat/GLAS receiver, and SLR data from the International Laser Ranging Service 
(ILRS) are needed. 


6.2 POD Products 

Two types of POD products will be generated: the Rapid Reference Orbit 
(RRO) and the operational POD. The former product will be generated within 12-24 
hours for primarily internal use of assessing the operational orbit and verification 
support for mission planning. The operational POD will be generated within 14 days, 
possibly within 3 days, after accounting for problems identified in RRO (e.g. GPS 
satellite problems) and problems reported by IGS. This product will be used in 
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generating the altimetry standard data products, particularly level IB and level 2 
surface elevation products. 


6.3 ICESat/GLAS Orbit and Attitude 

During the first 30-150 days after launch, the ICESat/GLAS spacecraft 
will be operated in a calibration orbit, with an 8-day repeat ground-track interval and 
94-degree inclination. At some point during this period to be determined by 
calibration results, the orbit will be transitioned to a neighboring mission orbit at the 
same inclination, with a 183-day repeating ground track. The ICESat/GLAS 
operational scenarios and orbit parameters are summarized in Table 6.1. 


Table 6.1 ICESat/GLAS Orbit Parameters 


Mission Phase 

Expected 

Duration 

(days) 

Mean 

Altitude 

(km) 

Inclina- 

tion 

(deg) 

Eccen- 

tricity 

Ground Track 
Repeat Cycle 

S/C Checkout 

30 

600 

94 

0.001 

No requirement 

Calibration/ 

Validation 

31-150 

600 

94 

0.0013 

8 days/183 days 

Polar Mapping 

151-1220 

600 

94 

0.0013 

183 days with 
25 and 8 day sub-cycles 


The ICESat/GLAS spacecraft will operate in two attitude modes 
depending on the angular distance between the orbit plane and the Sun (P' angle). As 
shown in Figure 1, for low-P' periods, such as that immediately following launch, the 
so-called "airplane-mode" is in use, with the solar panels perpendicular to the orbit 
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plane. When the P' angle exceeds 32 degrees, however a yaw maneuver places the 
satellite in the "sailboat-mode", with the axis of solar panels now in the orbit plane. 
While the two attitudes ensure that the solar arrays produce sufficient power year- 
round for bus and instrument operations, they introduce significantly different 
atmospheric drag effects due to the difference in cross-sectional area perpendicular to 
the velocity vector. 


6.4 POD Accuracy Assessment 

The predicted radial orbit errors based on recent gravity models (e.g., 
JGM-3 or EGM-96) are 19-36 cm. To reduce the effect of the geopotential model 
errors on ICESat/GLAS, which is the major source of orbit error for ICESat/GLAS 
POD, the gravity model improvement effort will be made through gravity tuning. 
Solar activity is predicted to peak shortly after launch, and decline significantly 
during the mission. The level of this activity correlates directly with the magnitude 
of atmospheric drag effects on the satellite. The combinations of high solar flux and 
low P' angle at the start of the mission poses special challenges for POD and gravity 
tuning. 

A previous simulation study [Rim et al., 1996] indicated that the 
ICESat/GLAS POD requirements could be met at 700-km altitude by either the 
gravity tuning or employing frequent estimation of empirical parameters, such as 
adjusting one-cycle-per-revolution parameters for every orbital revolution, within the 
context of a fully dynamic approach. This approach is referred to as a highly 
parameterized dynamic approach. Because the mission orbit altitude was lowered to 
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600-km, and the satellite design has been changed since this earlier study, a new in- 
depth simulation study [Rim et al., 1999] was conducted. It also indicates that even at 
600-km altitude with maximum solar activity, the 5-cm and 20-cm radial and 
horizontal ICESat/GLAS orbit detennination requirement can be met using this 
aforementioned gravity tuning and fully dynamical reduction strategy. Table 6.2 
summarizes the ICESat/GLAS orbit accuracy based on two geopotential models, pre- 
tune and post-tune models. The results are based on eight 1-day arcs with three 
different parameterizations. Those are (A) 1-rev Cd, 6-hour lcpr TN, (B) 1-rev Cd, 3- 
hour lcpr TN, and (C) 1-rev Cd, 1-rev lcpr TN, where 1-rev Cd indicates solving for 
drag coefficient for every orbital revolution, and lcpr TN means solving for one- 
cycle-per-revolution Transverse and Normal parameters. Note that even the case (C) 
could not meet the radial orbit detennination requirement using the pre-tune 
geopotential model. This indicates that gravity tuning is necessary to achieve the 
orbit determination requirement. A factor of three improvement in radial orbit 
accuracy was achieved for case (A), and a factor of two improvement occurred for 
case (C) by the post-tune gravity field. 


Table 6.2 ICESat/GLAS Orbit Errors (cm) 



Pre-Tune 

Post-Tune 

Case 

Data 

RMS Orbit Errors 

Data 

RMS Orbit Errors 


RMS 

R 

T 

N 

RMS 

R 

T 

N 

A 

5.0 

15.5 

35.2 

14.1 

1.9 

5.2 

11.2 

5.6 

B 

3.6 

10.3 

22.4 

11.2 

1.7 

3.6 

10.7 

5.4 

C 

2.3 

6.5 

12.2 

5.9 

1.6 

3.3 

10.1 

5.2 
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6.5 POD Processing Strategy 
6. 5. 1 Assumptions and Issues 

Several assumptions were made for the POD processing. We assume: 1) 
continued operation of IGS GPS network and the SLR network, 2) IGS GPS data is 
available in RINEX (Receiver Independent Exchange) format, 3) ICESat/GLAS GPS 
receiver has performance characteristics comparable to the flight TurboRogue, and 
ICESat/GLAS GPS data are available in RINEX format, and 4) most relevant IGS, 
SLR and ICESat/GLAS data are available within 24-36 hours. There are several 
issues for POD processing which include: 1) identification of problem GPS satellites, 
2) identification of problems with ground station data, 3) processing arc length, 4) 
accommodation for orbit maneuvers, and 5) problems associated with expected out- 
gassing during early mission phase. For a July 2001 launch and the early phases of 
the mission, orbit maneuvers are expected to occur as frequently as 5 days because of 
high level of solar activity [Demarest and Schutz, 1999]. These maneuvers will not 
be modeled, but the maneuver times will be utilized to reinitialize the orbit arc length. 


6. 5. 2 GPS Data Preprocessing 

The GPS data processing procedure consists of two major steps: data 
preprocessing and data reduction. The data preprocessing step includes data 
acquisition, correcting measurement time tags, generating double-differenced 
observables, and data editing. The GPS data preprocessing system is collectively 
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called TEX GAP (university of TEXas Gps Analysis Program) and implemented on 
the HP workstation. 

The International GPS Service for Geodynamics (IGS) provides GPS data 
collected from globally distributed GPS tracking sites, which include more than 200 
ground stations at present [IGS, 1998]. The daily IGS data files are archived in the 
IGS global data centers in the RINEX fonnat, and the data from selected ground 
station network will be downloaded to GSR’s data archive system. Also, the GPS 
data from the ICESat/GLAS GPS receiver will be provided by the ICESat Science 
Investigator Processing System. 

The GPS receiver time tag is in error due to the receiver clock error, and 
the time tag correction, t r , can be obtained by 


t r = p/C - p c /C + t s (6. 1) 

where C is the speed of light, p is the pseudorange measurement, p c is the computed 
range from GPS ephemerides and receiver position, and t s is the broadcast GPS 
satellite clock correction. 

Double-differencing eliminates common errors, such as the GPS satellite 
and receiver clock errors, including the Selective Availability (SA) effect. As 
described in Section 4.2.3, a double-differenced high-low observation consists of a 
ground station, two GPS satellites, and ICESat/GLAS satellite. A careful selection of 
double-differenced combination is required to avoid generating dependent data set. 

To eliminate the first-order ionospheric effects, the double-differenced 
carrier phase observables DD L1 at Li and DD L 2 at L 2 frequency are combined to form 
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where//./ = 1575.42 MHz and /12 = 1227.60 MHz. 

Data editing involves the detection and fixing of the cycle-slips of the 
carrier phase data, and the editing of data outliers. For editing outliers, a 3a editing 
criterion is applied to the double-differenced residual. Cycle-slips are detected by 
examining the differences between the consecutive data points in the double- 
differenced residuals and identifying discontinuity. The identified cycle-slips are 
fixed by using linear extrapolations. 

6. 5. 3 GPS Orbit Determination 

ICESat/GLAS POD requires precise GPS ephemerides, and there are two 
approaches to obtain the precise GPS ephemerides. The first approach is to solve the 
GPS orbit simultaneously with the ICESat/GLAS orbit, and the second approach is to 
fix the GPS ephemeris to an independent determination, such as the IGS solutions. 
For the first approach, standard models described in Table 6.3 will be used for the 
reference frame and gravitational perturbations for GPS. For the non-gravitational 
perturbations on GPS, the models described in Section 3.4.5 will be employed. It has 
been shown for the Topex POD case that adjusting GPS orbits usually resulted better 
Topex orbit solutions [Rim et al., 1995]. A simulation study [Rim et al., 2000b] 
indicates that fixing GPS orbits to high accuracy solutions would generate reasonably 
well-tuned gravity field, thereby, the POD accuracy requirement could be met with 
fixing GPS approach. As the accuracy of IGS solutions improved significantly 
[Kouba et al., 1998], fixing GPS ephemeris to IGS solutions would be a preferred 
approach for ICESat/GLAS POD. These two approaches will be evaluated using 
available tracking data during the pre-launch period, such as CHAMP and JASON, 
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and ICESat/GLAS tracking data during the vcrification/validation period. CHAMP 
POD accuracy was assessed when the GPS ephemeris is fixed to IGS solutions, such 
as the ultra-rapid, rapid, and final solutions [Rim et al., 2002a] 

6.5.4 Estimation Strategy 

The adopted estimation strategy for ICESat/GLAS POD is the dynamic 
approach with tuning of model parameters, especially the geopotential parameters. 
Simulation studies indicate that frequent estimation of empirical parameters is an 
effective way of reducing orbit errors. The solutions from the sequential filter with 
process-noise will be investigated as a validation tool for the highly parameterized 
dynamic solutions. Results of Davis [1996] and Rim et al. [2000a] show that both 
highly parameterized dynamic approach with gravity tuning and the reduced-dynamic 
approach yield comparable results in high fidelity simulations. This comparison will 
continue with the flight data. 


6.6 POD Plans 

This section describes planned POD activities during the pre-launch and 
the post-launch periods. 

6. 6. 1 Pre-Launch POD Activities 

During the pre-launch period, POD activities will be focused on the 
following areas: 1) selection of POD standards, 2) model improvement efforts, 3) 
preparation for operational POD, and 4) POD accuracy assessment. In this section, 
pre-launch POD activities in these areas are summarized. 



78 


6. 6. 1.1 Standards 

The standard models for the reference system, the force models and the 
measurement models to be used for the ICESat/GLAS POD are described in Table 
6.3. These standards are based on the International Earth Rotation Service (IERS) 
Conventions [ McCarthy , 1996], and the T/P standards [Tapley et ai, 1994]. These 
standards will be updated as the models improve, and “best” available models at 
launch will be selected as the initial standard models. 
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Table 6.3 Precision Orbit Determination Standards for ICESat/GLAS 

Model ICESat/GLAS Standard Reference 


Conventional inertial system 

Reference Frame 
ICRF 

IERS 

Precession 

1976 IAU 

IERS 

Nutation 

1980 IAU 

IERS 

Planetary ephemerides 

JPL DE-405 

Standish [1998] 

Polar Motion 

IERS 


UT1-TAI 

IERS 


Station Coordinates 

ITRF 


Plate motion 

Nuvel (NNR) 

IERS 

Reference ellipsoid 

a e = 6378136.3 m 

Wakker [1990] 

GM 

1/f = 298.257 

Force Models 
398600.4415 knf/s 2 

Ries et al. [1992a] 

Geopotential 

JGM-3 

Tapley et al. [1996] 


or EGM-96 

Lemoine et al. [1996] 


or TEG-4 

Tapley et al. [2001] 

C 21 , S 21 - mean values 

C 21 =-0.187xl0' 9 


C 21 , S 21 - rates 

5 2 i =+1.195xl0" 9 
C 21 = -1.3xl0" n /yr 

(see rotational deformation) 

Zonal rates 

5 2 i =+l.lxl0' u /yr 
epoch 1986.0 
J 2 = - 2.6x1 O' 1 Vyr 

Nerem et al. [1993] 

N body 

epoch 1986.0 
JPL DE-405 

Standish [1998] 

Indirect oblateness 

point mass Moon on Earth J 2 


Solid Earth tides 


IERS-Wahr [1981] 

Frequency independent 

k 2 = 0.3, k : s = 0.093 


Frequency dependent 

Wahr's theory 


Ocean tides 

CSR TOPEX 3.0 

Eanes and Bettadpur [1995] 

Rotational deformation 

AC21 = -1.3 xlO 9 (x p - x p ) 

Nerem et al. [1994] 

Relativity 

AS 2l =+1.3x10 -\y p -y p ) 

based on k 2 /ko = 0.319 
x p = 0".046, I P = 0".294 

Ip = 0".0033/yr 

y p = 0".0026/yr, epoch 1986.0 

all geocentric effects 

Ries et al. [1991] 

Solar radiation 

solar constant = 4.560x 10~ 6 



N hn 2 at 1 AU, conical shadow 
model for Earth and Moon 
R e = 6402 km, 

R, n = 1738 km, 
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R s = 696,000 km 


Atmospheric drag 

density temperature model 

Barlier et al. [1977] 


or MSIS90 

Hedin [1991] 


or NRLMSISE-00 
daily flux and 3 -hour 
constant k p , 

3-hour lag for k p , 

1-day lag for /10.7 , 

/io.7 average of previous 81 days 

Hedin et al. [1996] 

Earth radiation pressure 

Albedo and infrared second-degree 
zonal model, R e = 6378136.3 m 

Knocke et al. [1989] 

Satellite parameters 

ICESat/GLAS models 
Box-wing model 

Measurement Models 


Laser range 

Troposphere 

Marini & Murray [1973] 

IERS 

Relativity correction 

applied 

IERS 

Center of Mass/phase center 

ICESat/GLAS model 


GPS 

Troposphere 

MTT 

Herring [1992] 

Ionosphere 

dual frequency correction 


Center of Mass/phase center 

ICESat/GLAS model 


Relativity correction 

applied 


Site displacement 

Induced permanent tide 

IERS 


Geometric tides 

Frequency independent 

hi =0.6090, 
h = 0.0852, 
<5 = 0° 

IERS 

Frequency dependent 

K\ only 

IERS 

Ocean loading 

IERS 


Rotational deformation 

hi = 0.6090, h = 0.0852 with 
x p = 0".046, y p = 0",294 

x p = 0".0033/yr 

y p = 0".0026/yr, epoch 1986.0 

IERS 
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Figures 2 and 3 show the ground station network for ICESat/GLAS POD 
for GPS and SLR, respectively. Details of the adopted network may change prior to 
launch but will remain quite robust. Station coordinates will be adopted from the 
"best" available ITRF model, expected to be ITRF-99 or ITRF-2000. The ITRF 
model includes station velocities measured by space geodetic methods. 

6. 6. 1.2 Gravity Model Improvements 

The gravity model to be used in the immediate post-launch period will be 
"best" available at launch, such as JGM-3 [ Tapley et al., 1996], EGM-96 [Lemoine et 
al., 1996], or TEG-4 [ Tapley et ah, 2001]. As further gravity model improvements 
are made from other projects, such as GRACE, they will be incorporated for 
ICESat/GLAS POD. At this writing, further study is required for the selection of the 
at-launch gravity model. However, current state-of-the-art models are sufficiently 
close that geopotential tuning with ICESat/GLAS data should yield comparable POD 
performance which is largely unaffected by this initial selection. The “best” available 
ocean tide model at launch will be adopted as the standard ocean tide model for 
ICESat/GLAS POD. 

6. 6. 1. 3 Non-Gravitational Model Improvements 

Since the ICESat/GLAS launch coincides with the predicted solar 
maximum, the atmospheric drag perturbation will be the largest non-gravitational 
force acting on the satellite. Some drag-related models were evaluated for CHAMP 
POD, as part of drag model improvement efforts for reducing the effect of drag model 
errors on ICESat/GLAS POD [Rim et al., 2002b]. Those include the thermospheric 
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wind model, HWM93, NRLMSISE-00 [Hedin et al., 1996], and DTM-2000 
[Bruinsma and Thuillier, 2000]. Estimation strategies to minimize the effects of drag 
model errors on POD and gravity tuning will also be investigated. 

In order of decreasing magnitude, the remaining non-gravitational 
perturbations consist of solar radiation pressure, Earth radiation pressure, and on- 
board thermal emission. For POD, a ’box-wing’ model, described in Section 3.4.6, 
represents the spacecraft as a simple combination of a six-sided box and two attached 
panels, or 'wings’. This macro-model will use effective specular and diffuse 
reflectivity coefficients to compute the induced forces acting on each surface. The 
pre-flight values of these coefficients will be estimated during a tuning process, in 
which the forces computed with the macro-model are fit to those obtained using a 
separate micro-model [Webb, private communication, 2000], This latter model 
employs considerable detail that makes it impractical for use directly in POD. Once 
ICESat/GLAS is in orbit, the reflectivity coefficients will be adjusted during POD, 
using the GPS tracking data. 

The macro-model tuning effort will compute the radiation from various 
sources incident on the satellite's surfaces. By using a comprehensive thennal model, 
the propagation of this energy throughout the spacecraft will be calculated. The 
resulting temperature distribution will be evaluated to determine whether any on- 
board thermal gradients may induce net forces. Any such forces would then be 
modeled analytically during POD. 


The non-gravitational forces acting on each surface due to atmospheric 
drag, solar radiation pressure, Earth radiation pressure, and thennal emission are 
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computed individually and then summed to obtain the total non-gravitational force 
acting on the satellite. 

6. 6. 1. 4 Measurement Model Developments 

One of the sources of measurement model errors is the multipath effect. 
Colorado Center for Astrodynamics Research (CCAR) multipath study [ Axelrad et 
al., 1999] indicates that the multipath effect alone results in 1-2 cm radial orbit error, 
while this effect in the presence of other errors, such as drag and gravitational model 
errors, results in a few mm error. This study was based on a preliminary design 
location for the antennas and most of the multipath effect was caused by the solar 
arrays. It also indicates that the effect becomes even smaller with proper editing 
scheme, such as blocking certain regions. The capability of screening out GPS 
measurements from blocked regions was implemented in MSODP1. Strategies for 
detecting and mitigating the multipath effect on CHAMP POD were investigated 
[Yoon et al., 2002b], and similar approach will be adopted for ICESat/GLAS POD. 
The final spacecraft design has the GPS antennas positioned above the solar array and 
bus star cameras. In this location, there is no expected impingement above the 
ground plane so multipath will be mitigated. 

ICESat/GLAS satellite’s center of mass location with respect to a 
reference point on the spacecraft will be measured in the pre-launch period, and the 
location of the GPS antenna and the laser reflector will also be measured. GPS 
antenna phase center variations as a function of azimuth and elevation will be 
determined in pre-launch testing. Effect of GPS antenna phase center variation on 
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POD was investigated using CHAMP data [ Yoon et al., 2002a]. Expenditure of fuel 
and corresponding changes in center of mass location will be monitored during flight. 

6. 6. 1. 5 Preparation for Operational POD 

To generate the POD products operationally when large volumes of data 
are required, it is essential to make the POD processing as automatic as possible. The 
POD processing procedures will be examined end-to-end to identify/update the 
procedures for possible improvement and to minimize the human intervention, and 
computational and human resources will be allocated optimally for POD processing. 
The adopted operational POD processing procedures/scripts will be tested by 
processing upcoming satellites, such as JASON and CHAMP, during the pre-launch 
period for further improvement. 

6.6. 1.6 Software Comparison 

Since the POD products from different software will be compared for 
POD validation, it is important to compare different software packages in the pre- 
launch period to identify model differences and to quantify the level of agreement 
among different POD software systems, such as UT-CSR’s MSODP1, GSFC’s 
GEODYN, and JPL’s GOA II. This comparison becomes easier for the 
ICESat/GLAS POD due to the extensive POD software comparison activity between 
UT-CSR and GSFC for Topex POD [Ries, 1992b]. Also, Topex-GPS POD 
experiments between UT-CSR and JPL [ Bertiger et al., 1994] gave the opportunity 
for both groups to compare their software systems. This comparison will continue for 
the ICESat/GLAS POD models to ensure the validity of the POD verification by 
comparing with POD products from different software systems. 
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6. 6. 1. 7 POD Accuracy Assessment 

During the pre-launch period, simulation studies will continue to assess 
the POD accuracy. Comparison of highly parameterized dynamic approach and the 
reduced-dynamic approach will be continued. For the GPS orbit modeling, standard 
models for GPS orbit determination will be updated as the models progress, and the 
resulting orbit will be compared to the IGS solutions. Also, the effect of fixing GPS 
orbits to independently detennined ephemerides, such as IGS solutions, on the POD 
and the gravity tuning will be evaluated. 

6. 6. 2 Post-Launch POD Activities 

During the first 30-150 days after launch, which is the 
Calibration/Validation period, POD processing will tune the model parameters, 
including the gravity, and define adopted parameter set for processing the first 183- 
day cycle. During the 183 days of the Cycle 1, the POD processing will assess and 
possibly further improve or refine parameters, such as assess the gravity field from 
the gravity mission GRACE, if available, and adopt a new parameter set for the 
processing of Cycle 2 data. POD processing will continue assessment of POD quality 
after Cycle 1 , and new parameter adoptions should be minimized and timed to occur 
at cycle boundaries. 

6. 6. 2. 1 Verification/V alidation Period 

During the calibration/validation period, several important POD activities 
will be undertaken simultaneously. These include tuning model parameters, POD 
calibration/validation, evaluation of out-gassing effect, evaluation of estimation 
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strategies and GPS orbit modeling procedure, evaluation of multipath effect and 
construction of editing scheme. 

Some model parameters, such as geopotential parameters and the “box- 
wing” model parameters, will be tuned using the tracking data. About 30-40 days of 
GPS data will be processed for gravity tuning, and the arc length will be dictated by 
the maneuver spacing and the ability of POD to mitigate the effect of the non- 
gravitational model errors, especially the drag model errors, to certain level. The 
tuned gravity field will be determined by combining the pre-tune gravity coefficients 
and the solution covariance with the new information equations from the GPS 
tracking data. 

Internal and external POD calibration/validation activities are planned for 
POD quality assessment, and those are summarized in the following section. 

During the early phase of the mission, the satellite might experience 
significant out-gassing, and this poses serious challenges for POD. However, this 
effect will subside as time goes by, and every effort will be made to insure that this 
effect does not corrupt the parameter tuning process during this 
validation/verification period. 

Estimation strategies described in Section 6.5.4 will be evaluated, and the 
GPS orbit modeling procedures described in Section 6.5.3 will also be evaluated 
during this period. 

The multipath effect will be evaluated to characterize the extent of signal 
corruption due to diffraction and reflection using the flight data. Proper editing 
scheme will be developed if there is any evidence that such an editing reduces the 
multipath effect on POD. 
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6. 6. 2. 2 POD Product Validation 

To validate the accuracy of ICESat/GLAS POD products, several methods 
would be employed. For the internal evaluation of the orbit consistency, orbit 
overlap statistics will be analyzed. Also, the data fit RMS value is an effective 
indicator of orbit quality. Comparisons between the orbits from different software, 
such as MSODP1, GEODYN, and GIPSY-OASIS II (GOA II), would serve as a 
valuable tool to assess the orbit accuracy. Since the ICESat/GLAS will carry the 
laser reflector on board, the SLR data can be used as an independent data set to 
determine the ICESat/GLAS orbit. However, this approach assumes reasonably good 
tracking of the ICESat/GLAS orbit from the SLR stations. Data from the SLR 
network will also be used to directly evaluate the GPS-determined orbit. Data fits for 
high elevation SLR passes can be used to evaluate the orbit accuracy of the 
ICESat/GLAS. The laser altimeter data will be used to assess the validation, 
however, this assessment can be accomplished only if the calibration and verification 
of the instrument have been accomplished. Global crossovers from ICESat/GLAS 
will be used to validate the radial orbit accuracy in a relative sense. 

6. 6. 2. 3 POD Reprocessing 

To produce improved orbits, reprocessing of data will be performed as 
often as annually. As the solar activity is expected to decrease in the later mission 
period, the accuracy of the tuned model parameters will be improved, thereby the 
POD accuracy will be improved. Any improvement in the model parameters will be 
adopted for the reprocessing. 



6.7 Computational: CPU, Memory and Disk Storage 


Table 6.4 compares the computational requirements for processing a 
typical one-day arc from a 24-ground station network with 30 sec sampling time for 
both T/P and ICESat/GLAS. These results are based on MSODP1 implemented on 
the Cray J90 and the HP-735/125. 

Current computational plans are to use the HP-class workstation 
environment for preprocessing GPS data, including generation of double difference 
files. POD processing will be performed on a Cray J90, or equivalent. This 
processing on the Cray enables a more efficient resource sharing with other project, 
such as GRACE. 


Table 6.4 Computational Requirements for T/P and ICESat/GLAS POD using 


MSODP1: One-day Arcs with 24 

1 Ground Stations 

Platform 

Satellite 

CPU (min) 

Memory (Mw) 

Disk* (Mb) 

Cray J90 

T/P 

20 

2 

35 

ICESat/GLAS 

40 

2.5 

59 

HP-735 

T/P 

30 

2 

39 

ICESat/GLAS 

105 

2.5 

63 


This includes all the necessary files. 





























Figure 2. GPS Tracking Stations for ICESat/GLAS POD 



Figure 3. SLR Stations Tracking ICESat/GLAS (20 degree Elevation Masks) 
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Appendix A: ATBD Update for the Operational (“Final”) POD 
A.l ICESat Mission Outline 

The Ice, Cloud and land Elevation Satellite (ICESat) was launched on 13 
January 2003. The Geoscience Laser Altimeter System (GLAS) instrument onboard 
ICESat made its first laser elevation measurement of the Earth on 21 February 2003 
and its last on 1 1 October 2009. The three lasers employed by GLAS did not perfonn 
as long as expected, and following the failure of Laser 1 on 5 March 2003 the ICESat 
mission was modified to meet the requirement for capturing a multi-year time series 
of ice sheet elevations [ Schutz et al., 2005]. For the modified mission scenario, the 
spacecraft entered a 91 -day repeat science orbit (compared to a planned 183-day 
repeat) and the lasers were activated for about 33 days of this 91 -day repeat, two or 
three times per year. This campaign mode operation is summarized in Table A.l, and 
other significant parameters and events are listed in Table A. 2. ICESat laser 
campaigns are designated by a laser number (Li, L2 or L3), followed by a letter in 
the sequence of operation. Following campaign L2f, attempts to restart any of the 
lasers were not successful. The spacecraft was put through a series of engineering 
tests in early 2010. De-orbit maneuvers were carried out in June and July 2010. The 
spacecraft was “passivated” on 14 August and reentered the Earth’s atmosphere on 30 
August 2010 over the Barents Sea northeast of Norway. 
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Table A.l: ICESat Laser Operation Campaigns 


Campaign 

Year 

Day of year 

Calendar Dates 

Number of 
days (d) 

Repeat 
orbit (d) 

Repeat tracks 1 

Lla 

2003 

051-088 

20 Feb-29 Mar 

37 

8 

001-072 to 
006-023 

L2a 

2003 

268-277/ 

277-322 

25 Sep-4 Oct/ 
4 Oct- 18 Nov 

54 

8/ 

91 

028-088 to 029-100/ 
1098 to 0421 

L2b 

2004 

048-081 

17 Feb-21 Mar 

33 

91 

1284 to 0421 

L2c 

2004 

139-173 

18 May-21 Jun 

34 

91 

1283 to 0434 

L3a 

2004 

277-313 

3 Oct-8 Nov 

37 

91 

1273 to 0452 

L3b 

2005 

048-083 

17 Feb-24 Mar 

35 

91 

1258 to 0426 

L3c 

2005 

140-174 

20 May-23 Jun 

34 

91 

1275 to 0421 

L3d 

2005 

294-328 

21 Oct-24 Nov 

34 

91 

1282 to 0421 

L3e 

2006 

053-087 

22 Feb-28 Mar 

34 

91 

1283 to 0424 

L3f 

2006 

144-177 

24 May-26 Jun 

33 

91 

1283 to 0421 

L3g 

2006 

298-331 

25 Oct-27 Nov 

33 

91 

1283 to 0423 

L3h 

2007 

071-104 

12 Mar- 14 Apr 

33 

91 

1279 to 0426 

L3i 

2007 

275-309 

2 Oct-5 Nov 

34 

91 

1280 to 0421 

L3j 

2008 

048-081 

17 Feb-21 Mar 

33 

91 

1282 to 0422 

L3k 

2008 

278-293 

4 Oct- 19 Oct 

15 

91 

1283 to 0145 

L2d 

2008 

330-352 

25 Nov-17 Dec 

22 

91 

0096 to 0423 

L2e 

2009 

068-101 

9 Mar- 11 Apr 

33 

91 

1286 to 0424 

L2f 

2009 

273-284 

30 Sep-11 Oct 

11 

91 

1280 to 0084 


1 There are 119 tracks in the 8-day orbit and 1354 tracks in the 91 -day orbit. Cycle numbers 


are included for the 8-day repeat periods. 
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Table A.2: Significant ICESat Parameters and Events by Campaign 


Cam- 

paign 

Year 

Day 

of 

year 

S/C 

orient- 

ation 1 

Start 

Beta’ 

Angle 

(°) 

End 

Beta’ 

Angle 

(°) 

Start 

Laser 

Infrared 

Energy 

(mJ) 

End 

Laser 

Infrared 

Energy 

(mJ) 

Mean 
footprint 
major 
axis (m) 

Day of year - 
comments 

- 

2003 

013 

- 

- 

- 

- 

- 

- 

0 1 3 - launch 

Lla 

2003 

051-088 

-Y/+X 

-45 

-32 

72 

51 

149 

080 - yaw flip 

085 - safe hold, adjust 

temperature 

L2a 

2003 

268-277/ 

277-322 

+Y 

51 

69 

80 

55 

100 

277 - orbit change 
286 - laser temperature 
anomaly 

287, 302 - adjust 

temperature 

311 - GPS solar flare 

anomaly 

L2b 

2004 

048-081 

+Y 

54 

40 

57 

33 

90 


L2c 

2004 

139-173 

-X 

13 

-4 

33 

5 

88 

142-147 - adjust 
temperature 

L3a 

2004 

277-313 

-Y 

-48 

-58 

67 

62 

56 

293 - adjust 
temperature 

L3b 

2005 

048-083 

-Y 

-56 

-45 

68 

54 

80 

054 - suspected 
amplifier bar drop, 
begin footprint 
anomaly 2 
068 - suspected 
amplifier bar drop 

L3c 

2005 

140-174 

+X 

-20 

-4 

49 

44 

55 


L3d 

2005 

294-328 

+Y 

51 

63 

43 

39 

52 


L3e 

2006 

053-087 

+Y 

62 

48 

38 

30 

52 


L3f 

2006 

144-177 

-X 

20 

4 

30 

30 

51 

149 - Energy jump up 
2mJ 

L3g 

2006 

298-331 

-Y 

-44 

-54 

30 

24 

53 

310 -begin ITRF 2005 

L3h 

2007 

071-104 

-Y 

-60 

-47 

24 

21 

56 


L3i 

2007 

275-309 

+Y 

32 

46 

22 

20 

57 


L3j 

2008 

048-081 

+Y 

74 

62 

20 

16 

59 


L3k 

2008 

278-293 

+X 

-28 

-32 

18 

12 

52 

289 - Energy drop 4 
mJ 

L2d 

2008 

330-352 

-Y 

-45 

-53 

8 

4 

- 

343-344 - adjust 
temperature, energy up 
5 mJ 

L2e 

2009 

068-101 

-Y 

-71 

-59 

6 

2 

- 

094-095 - adjust 
temperature 

L2f 

2009 

273-284 

-X 

20 

25 

4 

2 

- 


- 

2010 

242 

- 

- 

- 

- 

- 

- 

242 - reentry 


1 The spacecraft is said to be in “Sailboat” mode for ±Y orientations and in “Airplane” mode 


for ±X orientations, where the direction indicates the solar panel orientation with respect to 
the spacecraft velocity using the GLAS coordinate frame. 

2 The footprint diameter during L3b changed from a mean of 54 m (day of year 048-053) to 
84 m (055-068). The reason for the larger footprint size during the latter part of the campaign 
is unknown, although a suspected amplifier bar dropout occurs near the event. 



94 


A.2 Gravitational Models 

It had been shown from pre-launch POD studies [Rim et al., 1996; Rim et al., 
1999] that the gravity model error was expected to be the dominant source of ICESat 
orbit errors. The predicted radial orbit errors at the ICESat orbit based on pre-launch 
gravity models, such as TEG-4 [ Tapley et al., 2000] and EGM-96 [Lemoine et al., 
1996], were 7-15 cm. As a consequence, gravity improvement using ICESat data 
(gravity tuning) was expected to be required. However, a gravity model from GRACE, 
GGM01C [ Tapley et al., 2004], was made available at ICESat launch, and the 
predicted radial orbit errors at the ICESat orbit was significantly reduced to 1 . 1 cm 
from this field. Since the gravity model error is no longer the dominant source of 
ICESat orbit error, the planned gravity tuning was not carried out. Note that 
GGM01C has been used for ICESat POD throughout the ICESat mission for the 
consistency of the POD products. Also note that CSR TOPEX 4.0 [Eanes and 
Bettadpur, 1995] was selected as the ocean tide model. 

A.3 Macro Model Development 

There has been a substantial effort to develop solar and Earth radiation 
pressure models for ICESat POD [ Webb et al., 2001; Rim et al., 2006]. ICESat 
macro-model [Webb, 2007] is the outcome for this endeavor. The model developed to 
compute solar and Earth radiation forces for ICESat POD consists of a six-sided box 
and two flat plates, or wings, representing the body of the satellite and its solar arrays, 
respectively. Illustrated in Fig. A.l (a), it is intended to capture the primary large 
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scale effects of incident radiation on satellite motion, and thus, has been designated 
the macro-model. This implementation requires knowledge of the specular and 
diffuse reflectivities, as well as the area, of each of the external surfaces of the box 
and both sides of each wing. Many types of geometric figures with varying materials, 
often with significantly different reflective properties, make up each face of the 
ICESat. Thus, to ensure that the macro-model adequately approximates the radiation 
pressure forces on ICESat, its parameters must effectively integrate the contributions 
from these disparate surfaces. 



Figure A.l. ICESat macro-model (a) and micro-model (b), with their common, 

body-fixed coordinate frame 

Adapting and extending the approach developed for the TOPEX/Poseidon 
mission [Antreasian and Rosborough, 1992; Marshall and Luthcke, 1992], the forces 
induced by solar, Earth albedo and Earth infrared radiation were simulated prior to 
the launch of ICESat, using a detailed model of the satellite, shown in Fig. A. 1 (b). 
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Ball Aerospace originally developed this micro-model with the Thennal Synthesizer 
System (TSS) for its thennal analyses of the satellite bus. It consists of 950 surfaces, 
including flat plates, cones and cylinders. Many of these are further subdivided, 
yielding 1124 nodes that can receive external radiation. As designed, the TSS 
software computes the incident and absorbed heat rates at each of these nodes using a 
sophisticated Monte Carlo Ray Tracing method. In consultation with the Center for 
Space Research, at the University of Texas at Austin, the TSS vendor, Space Design, 
Inc., implemented a series of modifications to employ a similar technique in the 
detennination of radiation pressure forces [ Webb et ah, 2001]. 

Using this enhanced software, the radiation forces acting at each node in the 
micro-model were computed at discrete points around the orbit. Once summed to 
obtain the net force components in the satellite body-fixed frame, they were rotated to 
an orbit coordinate system and expressed in radial, transverse, and normal (RTN) 
components. Generated at /?' angles every 5° between -90° and +90°, these single- 
revolution force histories collectively constituted a set of “truth” observations 
spanning the various orbit-Sun geometries and satellite orientations expected 
throughout the mission. These data were then fit using a least-squares (LSEI) method 
[ Hanson and Haskell, 1982] to find the best set of macro-model parameters for 
ICESat POD. This particular approach incorporated linear equality and inequality 
constraints to avoid physically unrealistic estimates. To ensure that no errors would 
inadvertently be introduced, the orbit, attitude, coordinate-transfonnation and Sun- 
position data were read from files output during the micro-model simulation. The 
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results, which were adopted for POD use during the mission, are shown in Table A. 3. 
Note that the coordinate set (X, Y, Z) in Table A. 3 is the Satellite Body-Fixed 
Coordinate System (SBCS), where X-axis is in the zenith direction and Y-axis is 
along the solar panel axis. The origin of SBCS is on the center longitudinal axis and 
8" above separation plane. Note also that YSA stands for Solar Array axis of rotation, 
which is in Y-direction in SBCS. 

Table A.3. ICESat Macro-Model Parameters 


Macro-Model 

Face 

Surface Area 
(m 2 ) 

Specular 

Reflectivity 

Diffuse 

Reflectivity 

+X 

3.82 

0.491 

0.000 

-X 

3.82 

0.951 

0.000 

+Y 

5.21 

0.426 

0.000 

-Y 

5.21 

0.413 

0.000 

+Z 

2.73 

0.345 

0.170 

-z 

2.73 

0.736 

0.000 

+YSA (front) 

4.21 

0.258 

0.000 

+YSA (back) 

4.21 

0.557 

0.000 

-YSA (front) 

4.21 

0.258 

0.000 

-YSA (back) 

4.21 

0.557 

0.000 


A.4 GPS antennae and Laser Retro-reflector Array (LRA) location measurement 

In the pre-launch period, Ball Aerospace & Technologies Corporation, which 
built the ICESat spacecraft Bus, measured the location of GPS antennae and LRA 
[Iacometti, 2002]. In Table A.4 the measured locations of GPS antennae are listed. 
Note that the measured GPS antenna location is the center of choke ring outboard 
surface. The computed LRA phase center location is also given in the Table A.4. 
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Note that additional 4.5 cm LRA correction should be applied when SLR data is 
processed to model LRA phase center location [Ries, 2003]. 


Table A.4. Pre-launch Measurements for GPS Antennae and LRA in Spacecraft 

Body-Fixed Coordinate System 



X 

Y 

Z 

FM- 1 Antenna Reference Point 

1.313 

0.189 

0.586 

FM-2 Antenna Reference Point 

1.313 

-0.189 

0.586 

LRA 

-0.99247 

0.0 

1.273 


All units m 


Satellite mass and the location of the Center of Mass (CoM) are changing after each 
maneuver, and Table A. 5 lists the orbit maintenance maneuvers and the mass and the 
location of CoM during ICESat campaigns. 


A.5 Estimated Parameters 

Estimated parameters in orbit determination process include ICESat state at 
the arc epoch, drag scaling parameters (Cd) for every orbital revolution, sinusoidal 
along-track (AT) and cross-track (CT) forces with a period equal to the orbital period 
(i.e., one cycle per revolution or 1-cpr) for every orbital revolution, double- 
differenced ambiguity parameters, piecewise constant zenith delay parameters for 
every 2.5-hours for ground stations, and the radial-component of center-of-mass 
offset correction parameter for the operating ICESat GPS antenna. 
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Table A.5. Orbit Maintenance Maneuvers, Center of Mass Location, and Satellite Mass 


Camp- 

Maneu- 

Time 

AV 

Aa* 

Center of Mass (m) 

Mass 

aign 

uver # 

(mm/dd/yy hh:mm:sec) 

(m/s) 

(m) 

X 

Y 

z 

(kg) 


- 

- 

- 

- 

0.052 

0.004 

0.984 

951.827 


3 

02/27/03 17:34:45.2682 

0.043901 

80.9249 

0.052 

0.004 

0.985 

951.807 

Lla 

4 

03/05/03 17:47:33.2300 

0.064421 

118.7504 

0.052 

0.004 

0.985 

951.778 


5 

03/14/03 15:17:22.1124 

0.071190 

131.2296 

0.052 

0.004 

0.985 

951.746 


6 

03/22/03 23:53:55.3704 

0.047863 

88.2284 

0.052 

0.004 

0.985 

951.724 


- 

- 

- 

- 

0.052 

0.004 

0.987 

948.567 


32 

09/26/03 13:04:59.2871 

0.041459 

-76.4243 

0.052 

0.004 

0.987 

948.595 


33 

10/04/03 12:25:34.3384 

0.337591 

-622.3018 

0.052 

0.004 

0.987 

948.447 


34 

10/04/03 13:13:54.5188 

0.337591 

-622.3018 

0.052 

0.004 

0.987 

948.299 


35 

10/07/03 03:39:09.7195 

0.020607 

-37.9862 

0.052 

0.004 

0.988 

948.279 

L2a 

36 

10/19/03 01:34:50.7209 

0.051738 

95.3720 

0.052 

0.004 

0.988 

948.240 


37 

10/25/03 21:04:55.1758 

0.047051 

86.7315 

0.052 

0.004 

0.988 

948.299 


38 

11/01/03 03:41:55.5413 

0.103132 

190.1093 

0.052 

0.004 

0.988 

948.295 


39 

11/04/03 00:09:35.5138 

0.028300 

-52.1671 

0.052 

0.004 

0.988 

948.220 


40 

11/11/03 02:06:22.3126 

0.048364 

89.1530 

0.052 

0.004 

0.988 

948.220 


- 

- 

- 

- 

0.052 

0.004 

0.988 

946.804 

L2b 

51 

02/28/04 23:54:15.3933 

0.043696 

80.5466 

0.052 

0.004 

0.989 

946.739 


52 

03/12/04 20:23:05.6900 

0.040266 

74.2252 

0.052 

0.004 

0.988 

946.949 


- 

- 

- 

- 

0.052 

0.004 

0.989 

946.633 

L2c 

59 

05/28/04 09:34:17.3027 

0.037728 

69.5455 

0.052 

0.004 

0.989 

946.586 


60 

06/10/04 06:03:11.2784 

0.037603 

69.3159 

0.052 

0.004 

0.988 

946.692 


- 

- 

- 

- 

0.052 

0.004 

0.989 

945.335 

L3a 

74 

10/12/04 05:37:53.5177 

0.032631 

60.1504 

0.052 

0.004 

0.990 

945.307 


75 

10/27/04 20:57:43.0406 

0.038153 

70.3305 

0.052 

0.004 

0.989 

945.398 


- 

- 

- 

- 

0.052 

0.004 

0.990 

944.688 

L3b 

88 

02/21/05 07:30:17.4819 

0.049142 

90.5856 

0.052 

0.004 

0.990 

944.668 


89 

03/14/05 23:48:47.4269 

0.041237 

76.0145 

0.052 

0.004 

0.990 

944.656 


- 

- 

- 

- 

0.052 

0.004 

0.990 

943.805 

L3c 

96 

06/05/05 13:10:18.6167 

0.041758 

76.9753 

0.052 

0.004 

0.991 

943.795 


97 

06/20/05 10:47:00.5638 

0.027876 

51.3854 

0.052 

0.004 

0.990 

943.916 


- 

- 

- 

- 

0.053 

0.004 

0.992 

941.544 

L3d 

Ill 

10/24/05 07:45:57.6404 

0.026643 

49.1131 

0.053 

0.004 

0.992 

941.911 


112 

11/12/05 18:54:55.8686 

0.022183 

40.8917 

0.053 

0.004 

0.992 

941.977 



- 


- 

0.053 

0.004 

0.993 

941.304 

L3e 

120 

03/12/06 03:50:43.5721 

0.027713 

51.0846 

0.053 

0.004 

0.992 

941.286 






0.053 

0.004 

0.993 

940.718 

L3f 

125 

06/07/06 21:14:28.2302 

0.025539 

47.0767 

0.053 

0.004 

0.993 

940.749 


- 

- 

- 

- 

0.053 

0.004 

0.994 

939.549 

L3g 

135 

11/17/06 19:25:04.9200 

0.023362 

43.0647 

0.053 

0.004 

0.994 

939.602 


- 

- 

- 

- 

0.053 

0.004 

0.994 

938.769 

L3h 

143 

03/23/07 15:28:55.9619 

0.013268 

24.4579 

0.053 

0.004 

0.994 

938.915 


144 

04/08/07 06:48:46.7778 

0.021437 

39.5160 

0.053 

0.004 

0.994 

938.698 

L3i 

no maneuvers 

0.053 

0.004 

0.995 

937.674 





- 

0.053 

0.004 

0.996 

936.567 

L3j 

166 

03/05/08 04:34:33.7579 

0.017996 

33.1731 

0.053 

0.004 

0.996 

936.547 

L3k 

no maneuvers 

0.053 

0.004 

0.996 

936.179 


- 

- 

- 

- 

0.053 

0.004 

0.998 

933.910 

L2d 

181 

12/14/08 07:44:48.9590 

0.016426 

30.2791 

0.053 

0.004 

0.997 

935.022 





- 

0.053 

0.004 

0.998 

934.104 

L2e 

188 

03/24/09 12:54:57.8440 

0.013900 

25.6219 

0.053 

0.004 

0.998 

934.361 





- 

0.053 

0.004 

0.999 

932.724 

L2f 

198 

10/054)9 21:29:10.0982 

0.017056 

31.4396 

0.053 

0.004 

0.999 

933.012 


* change in the semi-major axis 
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A.6 POD Processing Strategy 

It was suggested in the POD ATBD (version 2.2) that the reduced dynamic 
solutions will be investigated as a validation tool for the highly parameterized 
dynamic solutions. However, no reduced dynamic solutions were generated and 
compared with the dynamic solutions during ICESat mission. Results of Davis [1996] 
and Rim et al. [2000a] show that both highly parameterized dynamic approach with 
gravity tuning and the reduced-dynamic approach yield comparable results in high 
fidelity simulations. With the GGM01C field and favorable solar activities throughout 
ICESat mission, except in 2003, it is expected that the two approaches would 
generate comparable results. 

The GPS orbits were fixed to IGS orbits, and station coordinates were fixed to 
ITRF2000 solutions [ Altamimi et al., 2002] up to L3f campaign, and to ITRF2005 
solutions [Altamimi et al., 2007] starting with the L3g campaign. Note that the IGS 
switched from ITRF2000 to ITRF2005 during the L3g campaign (Nov 05, 2006) to 
generate the GPS orbits. ICESat attitude was modeled by PAD solutions, and the on- 
board solar array orientation infonnation was used. 

A.7 POD Accuracy Assessment 

A common method for evaluating the POD precision is to compare 
ephemerides that are adjacent in time and overlap. For ICESat POD, a 30-hour arc is 
processed, where the middle 24-hour portion is the daily POD product and the 
additional 6-hours overlaps with the independently detennined adjacent arcs before 
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and after the 24-hour product. Orbit comparison in the overlapping region provides a 
measure of precision. Such a measure is somewhat representative of internal precision, 
but it also is indicative of orbit accuracy, although it usually provides an optimistic 
estimate. Table A. 6 summarizes the mean double-differenced RMS for both Rapid 
and Final POD solutions. The double-differenced RMS for all solutions is about 1 cm. 


Table A.6. Mean DD-RMS 


Campaign 

Rapid POD 

Final POD 

Lla 

1.03 

1.01 

L2a 

1.07 

1.02 

L2b 

1.04 

1.01 

L2c 

1.06 

1.04 

L3a 

1.00 

1.00 

L3b 

1.01 

1.00 

L3c 

1.02 

1.01 

L3d 

1.03 

1.02 

L3e 

1.01 

1.01 

L3f 

1.07 

1.07 

L3g 

1.02 

1.02 

L3h 

1.04 

1.03 

L3i 

1.04 

1.03 

L3j 

1.01 

1.01 

L3k 

1.02 

1.01 

L2d 

1.05 

1.05 

L2e 

1.06 

1.05 

L2f 

1.10 

1.10 

Mean 

1.04 

1.03 


All units cm 


Table A. 7 shows the mean orbit overlap statistics for each campaign. The 
mean of the mean radial orbit overlap for all campaigns is less than 7 mm for both 
Rapid and Final POD solutions, and the mean of 3D RSS (3 -Dimensional Root Sum 









102 


Square) is less than 1 .4 cm. Note that there is slight improvement in the DD-RMS and 


orbit overlaps in the Final solutions comparing with the Rapid solutions. 


Table A.7. Mean Overlap Statistics 


Campaign 

Rapid 

IPOD 

Final POD 

R 

T 

N 

3D 

RSS 

R 

T 

N 

3D 

RSS 

Lla 

0.70 

1.16 

0.73 

1.56 

0.63 

1.03 

0.68 

1.42 

L2a 

0.82 

1.26 

0.75 

1.72 

0.79 

1.36 

0.67 

1.77 

L2b 

0.80 

1.27 

0.80 

1.74 

0.71 

1.01 

0.69 

1.46 

L2c 



0.72 

1.50 



0.70 

1.41 

L3a 



0.69 

1.47 



0.69 

1.41 

L3b 



0.81 

1.42 



0.76 

1.32 

L3c 



0.49 

1.32 



0.48 

1.32 

L3d 

0.60 

0.96 

0.56 

1.28 

0.59 

0.96 

0.57 

1.29 

L3e 

0.59 

0.93 

0.54 

1.24 

0.59 

0.89 

0.54 

1.21 

L3f 

0.67 

1.16 

0.77 

1.57 

0.66 

1.11 

0.78 

1.53 

L3g 

0.52 

0.82 

0.60 

1.16 

0.53 

0.85 

0.60 

1.19 

L3h 

0.56 

0.88 

0.69 

1.28 

0.55 

0.88 

0.73 

1.29 

L3i 

0.66 

0.93 

0.55 

1.29 

0.64 

0.90 

0.54 

1.24 

L3i 

0.51 

0.87 

0.61 

1.20 

0.51 

0.87 

0.62 

1.22 

L3k 

0.52 

0.75 

0.44 

1.04 

0.51 

0.73 

0.45 

1.03 

L2d 

0.64 

0.97 

0.55 

1.31 

0.61 

0.95 

0.56 

1.28 

L2e 

0.66 

1.03 

0.67 

1.41 

0.63 

0.97 

0.64 

1.35 

L2f 

0.57 

1.10 

0.54 

1.38 

0.58 

1.10 

0.51 

1.36 

Mean 

0.64 

1.01 

0.64 

1.38 

0.62 

0.97 

0.62 

1.34 


All units cm 


The ground-based laser ranging measurements provide an independent tool to 
assess the accuracy of ICESat POD. By withholding the SLR data from the POD 
solution, range residuals can be formed using the adopted tracking station coordinates 
and the GPS -determined ICESat ephemeris. Ten SLR stations participate in the 
ranging to ICESat: Zimmerwald, McDonald Observatory, Yarragadee, Greenbelt, 


Monument Peak, Haleakala, Graz, Herstmonceux, Arequipa, and Hartebeesthoek. To 
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assure the safety of the detector in the GLAS instrument, a 70-degree maximum 
elevation pointing restriction has been imposed on those ground stations participating 
in the ranging. This prevents ground based lasers from sending laser radiation into the 
GLAS telescope, which could potentially damage the laser detector used by GLAS. 

SLR residuals reflect not only the radial component of orbit errors, but also 
the horizontal component. Usually, the high elevation residuals represent more of the 
radial component of the orbit accuracy, but unfortunately, there was no tracking over 
70-degree due to the restriction mentioned above. To quantify the contribution of the 
radial orbit errors to the SLR residuals, a “radial” RMS was calculated based on the 
tracking geometry. Also, high elevation residuals are computed using SLR data above 
60 degree elevation. Note that ITRF-2005 SLR station coordinates were used for this 
analysis. 

Table A. 8 summarizes the SLR residual statistics for the Final POD. Note that 
the tracking from the participating stations was increased substantially as the ICESat 
mission progresses, and average of 6.6 passes were tracked per day for L3f through 
L3i campaigns. The overall RMS was 1.90 cm. The “radial” RMS was 1.19 cm, and 
the high elevation (above 60 degree) RMS was 1.65 cm. 

Table A. 9 summarizes the SLR residual statistics by the tracking stations for 
the Final POD. Zimmerwald tracks with two laser frequencies, one for infrared (I), 
and the other for violet (V). Yarragadee has the best tracking with 555 passes total. 
Hartebeesthoek generated the highest RMS of about 3.7 cm. Haleakala and 
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Yarragadee performed the best with RMS of about 1.6 cm, and the “radial” RMS of 
about 1 cm. Among all stations which have more than 100 passes of tracking, 
Yarragadee perfonned the best with 1.58 cm RMS and 0.97 cm “radial” RMS. 
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Tal 

ble A.8. 1 

ICESat SLR Residuals for Final POD 






Range 


RMS 

Range 

Campaign 

Pass # 

Data # 

RMS 

Bias 

Radial 

(data #) 
>60° 

Bias 

>60° 

Lla 

22 

330 

1.43 

1.01 

0.87 

n/a 

n/a 

L2a 

8 

272 

1.78 

1.03 

1.35 

1.82(10) 

0.53 

L2b 

3 

37 

1.76 

0.22 

0.95 

n/a 

n/a 

L2c 

23 

250 

1.62 

0.96 

1.03 

1.32(2) 

n/a 

L3a 

6 

48 

1.29 

0.78 

0.80 

n/a 

n/a 

L3b 

37 

999 

1.68 

0.97 

0.99 

1.50(13) 

0.76 

L3c 

68 

1860 

2.02 

0.90 

1.26 

1.52(51) 

1.16 

L3d 

191 

6278 

1.84 

1.06 

1.10 

1.56(142) 

1.43 

L3e 

85 

3146 

1.74 

1.20 

1.08 

1.60(65) 

1.82 

L3f 

216 

7741 

2.07 

1.15 

1.32 

1.75(311) 

1.52 

L3g 

241 

8803 

2.27 

1.08 

1.43 

1.88(346) 

1.52 

L3h 

222 

8561 

1.69 

0.89 

1.06 

1.45(304) 

1.38 

L3i 

235 

8270 

1.83 

1.01 

1.13 

1.42(267) 

1.09 

L3j 

115 

5103 

1.77 

1.05 

1.12 

1.57(225) 

1.37 

L3k 

99 

3813 

1.38 

1.01 

0.89 

1.48(127) 

1.44 

L2d 

86 

3308 

2.09 

1.21 

1.31 

2.09(113) 

1.71 

L2e 

178 

7113 

2.00 

1.16 

1.22 

1.58(279) 

1.40 

L2f 

53 

2290 

1.72 

1.19 

1.07 

1.94(54) 

1.66 

Total 

1888 

68222 

1.90 

1.09 

1.19 

1.65(2309) 

1.46 


All units cm 


Table A.9. ICESat SLR Residuals by Stations for Final POD 


Station 

Pass 

# 

Data # 

RMS 

Range 

Bias 

Radial 

RMS 
(data #) 
>60° 

Range 

Bias 

>60° 

Zimmerwald-I 

205 

5612 

2.33 

1.14 

1.47 

1.97(187) 

1.19 

Zimmerwald-V 

253 

8181 

2.01 

0.89 

1.25 

1.47(261) 

1.14 

McDonald 

86 

1147 

2.08 

1.60 

1.23 

n/a 

n/a 

McDonald2 

17 

316 

2.98 

2.15 

1.70 

n/a 

n/a 

Yarragadee 

555 

24847 

1.58 

1.07 

0.97 

1.62(697) 

1.49 

Greenbelt 

134 

5629 

1.88 

1.21 

1.08 

1.76(111) 

1.27 

Monument Peak 

144 

5072 

1.76 

1.20 

1.15 

1.71(213) 

1.14 

Haleakala 

27 

916 

1.58 

1.23 

1.02 

2.64( 6) 

0.00 

Graz 

227 

9186 

1.95 

0.88 

1.14 

1.37(257) 

0.94 

Herstmonceux 

171 

5729 

2.02 

1.05 

1.45 

1.68(523) 

1.38 

Arequipa 

23 

265 

1.68 

1.15 

1.14 

2.48( 9) 

2.58 

Hartebeesthoek 

46 

1322 

3.70 

0.98 

2.32 

2.04( 45) 

1.60 

Total 

1888 

68222 

1.90 

1.09 

1.19 

1.65(2309) 

1.46 


All units cm 
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Appendix B: 2011 POD Reprocessing 

In 2011, GPS data for all the campaign periods during ICESat mission was 
reprocessed to ensure consistency within ICESat POD products among each 
campaign period and to reflect POD model advancements since the Final POD 
standard had been selected. Table B.l summarizes the changes made for this 
reprocessing from Final POD. 


Table B.l. Changes in 2011 POD Reprocessing 


• POD hardware platform and OS: 

o HP Workstation => Workstation with Intel CPU and Linux OS 

• POD Fortran compiler: 

o HP Fortran Compiler => Intel Fortran Compiler 

• MSODP version: 

o MSODP 2003. 1 => MSODP 2009. 1 (EPHEVL and ROTATE from 
2003.1) 

• Reference frame: 

o IGS GPS ephemeris products: IGS final => IGS reprol (LI ~ L3i) 
o GPS ground station network: 42 common, 10 deleted, 3 added 
o GPS ground station location: IGS00 => IGS05 (LI ~ L3f), ITRF2005 
=> IGS05 (L3g ~ L2f) 

o EOPDAT : EOPDATC to EOPDATC05 (LI ~ L3h) 

• Gravitational models: 

o Geopotential field: GGM01C => GGM03C 

o Ocean tide model: CSR TOPEX 4.0 => FES2004 

o Solid Earth tide model: IERS 1996 => IERS 2003 

• Observation models: 

o Ground antenna phase center offset: igs_0 1 .pcv => igs05 .atx 

o GPS transmitter phase center variation: to igs05.atx 

o Updated GPS yaw table: for LI ~ L2d 

• Estimated parameters: 

o Cross-track COM correction parameter added 
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B.l POD Environment 

Final POD processing was perfonned using an HP workstation. Reprocessing 
was performed on a workstation with Intel CPU and Linux OS. The primary software 
for POD used in Final POD was MSODP version 2003.1 compiled with HP 
FORTRAN compiler. In reprocessing, MSODP version 2009.1 compiled with Intel 
FORTRAN compiler was used. 

B.2 Reference Frame 

IGS reprocessed all the data prior to 2008 using IGS05 \Gendt, 2006], and 
generated “reprol” products [ Ferland , 2010], and the reprocessing used these 
products for LI a through L3i campaigns. 

For Final POD, ground station network coordinates information was taken 
from the file IGS00 [Weber, 2001] for Lla through L3f campaigns and from the file 
ITRF2005 for L3g through L2f campaigns, respectively. In reprocessing, that 
information for all the campaigns was taken from the file IGS05. 

GPS ground station network was updated for reprocessing. Ten stations were 
deleted and three stations were added to the station network that was used for Final 
POD. There were 42 common stations between the old and the new station network. 

Concerning the Earth orientation file EOPDAT, for campaigns from Lla to 
L3h, the file EOPDAT C was used in Final POD. In reprocessing for these 
campaigns, EOPDAT C 05, which is consistent with ITRF 2005 reference frame, 
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was used. For campaigns from L3i to L2f, the file EOPDATC05 was used in Final 
POD as well as in reprocessing. 

B.3 Gravitational Models 

Gravity field was updated from GGM01C to GGM03C [Tapley et ai, 2007] 
and ocean tide model was updated from CSR TOPEX4.0 to FES2004 [ Lyard et al., 
2006] for the reprocessing. Solid Earth tide model was based on IERS 1996 
[ McCarthy (eds.), 1996] standard for Final POD, and this model was updated using 
IERS 2003 [McCarthy and Petit (eds.), 2004] standard for reprocessing. This change 
was introduced by the MSODP update from version 2003.1 to version 2009. 1 . 

B.4 Observation Models 

The ground antenna phase center offset information was taken from the file 
“igs Ol.pcv” for Final POD. In reprocessing, that infonnation was taken from the file 
“igs05.atx”. GPS transmitter phase center variation was not modeled in Final POD, 
but it was modeled based on “igs05.atx” file for reprocessing. 

B.5 Estimated Parameters 

In Final POD, the radial component of center-of-mass offset correction 
parameter for the operating ICESat GPS antenna was estimated once per each 
estimation arc. In reprocessing, the cross-track component as well as the radial 
component was estimated once per each estimation arc. 
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B.6 Reprocessed POD Accuracy Assessment 

Table B.2 summarizes the mean orbit differences between the Final POD and 
the reprocessed POD. Note that there was a period where there was no GPS tracking 
data due to the orbit event at the end of LI a campaign, and this period was excluded 
for the orbit comparison. Mean radial differences were about 6 mm, and the mean 
3D-RSS differences were about 1.4 cm. 


Table B.2. Orbit Difference between Final POD and the reprocessed POD 


Campaign 

R 

T 

N 

3D-RSS 

Lla 

0.59 

1.26 

0.70 

1.56 

L2a 

0.49 

1.00 

0.64 

1.29 

L2b 

0.55 

1.19 

0.60 

1.45 

L2c 

0.60 

1.24 

0.71 

1.56 

L3a 

0.50 

1.07 

0.62 

1.33 

L3b 

0.61 

1.20 

0.63 

1.49 

L3c 

0.57 

1.25 

0.67 

1.53 

L3d 

0.59 

1.27 

0.65 

1.55 

L3e 

0.58 

1.19 

0.64 

1.47 

L3f 

0.74 

1.43 

0.71 

1.76 

L3g 

0.50 

0.94 

0.64 

1.24 

L3h 

0.57 

1.07 

0.68 

1.39 

L3i 

0.56 

1.22 

0.66 

1.50 

L3j 

0.63 

1.18 

0.63 

1.48 

L3k 

0.50 

1.01 

0.66 

1.31 

L2d 

0.52 

1.03 

0.55 

1.28 

L2e 

0.58 

1.11 

0.56 

1.38 

L2f 

0.49 

1.15 

0.60 

1.39 

Mean 

0.57 

1.16 

0.64 

1.44 


All units cm 


Table B.3 compares the mean double-differenced RMS between the Final 
POD and the reprocessed POD. There was slight improvement in DD RMS for 
reprocessed POD in sub-millimeter level. Table B.4 compares the mean orbit 






110 


overlaps between the Final POD and the reprocessed POD. There was similar 
improvement in the mean orbit overlaps in sub-millimeter level for the reprocessed 
POD. 


Table B.3. Mean DD-RMS 


Campaign 

Final POD 

2011 Reprocessing 

Lla 

1.01 

1.00 

L2a 

1.02 

1.01 

L2b 

1.01 

0.99 

L2c 

1.04 

1.03 

L3a 

1.00 

0.98 

L3b 

1.00 

0.98 

L3c 

1.01 

1.00 

L3d 

1.02 

1.01 

L3e 

1.01 

0.99 

L3f 

1.07 

1.05 

L3g 

1.02 

0.99 

L3h 

1.03 

1.01 

L3i 

1.03 

1.00 

L3j 

1.01 

0.97 

L3k 

1.01 

0.97 

L2d 

1.05 

1.01 

L2e 

1.05 

1.02 

L2f 

1.10 

1.04 

Mean 

1.03 

1.00 


All units cm 


Table B.5 summarizes the SLR residual statistics for the reprocessed POD. 
The overall RMS was 1.75 cm, which is 1.5 mm smaller than the overall RMS for 
Final POD. The “radial” RMS was 1.09 cm, 1 mm improvement over the Final POD 
case, and the high elevation (above 60 degree) RMS was 1.58 cm, which is 0.7 mm 


smaller than the Final POD case. 
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Table B.6 summarizes the SLR residual statistics by the tracking stations for 
the reprocessed POD. The overall RMS for the reprocessed POD reduced for all 
stations, except Haleakala and Hartebeesthoek. Among all stations, Yarragadee 
performed the best with 1.47 cm RMS and 0.91 cm “radial” RMS. 


Table B.4. Mean Overlap Statistics 


Campaign 

Final POD 

20 1 1 Reprocessing 

R 

T 

N 

3D 

RSS 

R 

T 

N 

3D 

RSS 

Lla 

0.63 

1.03 

0.68 

1.42 

mm 


0.56 

1.22 

L2a 

0.79 

1.36 

0.67 

1.77 

mm, 


0.57 

1.62 

L2b 

0.71 

1.01 

0.69 

1.46 

mm 


0.48 

1.17 

L2c 

mm. 

1.01 

0.70 

1.41 

mm 

1.09 

mm. 

1.49 

L3a 

mm. 

1.02 

0.69 

1.41 

mm 

0.83 

imm 

1.18 

L3b 

0.54 

0.89 

0.76 

1.32 

0.61 

0.87 

0.65 

1.27 

L3c 

0.74 

0.94 

0.48 

1.32 

0.71 

0.94 

0.40 

1.28 

L3d 

0.59 

0.96 

0.57 

1.29 

0.55 

0.92 

0.49 

1.20 

L3e 

0.59 

0.89 

0.54 

1.21 

0.51 

0.75 

0.41 

1.02 

L3f 

0.66 

1.11 

0.78 

1.53 

0.75 

1.22 

0.66 

1.61 

L3g 

0.53 

0.85 

0.60 

1.19 

0.53 

0.82 

0.58 

1.15 

L3h 

0.55 

0.88 

0.73 

1.29 

0.59 

0.92 

0.67 

1.31 

L3i 

0.64 

0.90 

0.54 

1.24 

0.47 

0.76 

0.50 

1.05 

L3j 

0.51 

0.87 

0.62 

1.22 

0.46 

0.71 

0.53 

1.02 

L3k 

0.51 

0.73 

0.45 

1.03 

0.56 

0.78 

0.37 

1.05 

L2d 

0.61 

0.95 

0.56 

1.28 

0.50 

0.85 

0.50 

1.13 

L2e 

0.63 

0.97 

0.64 

1.35 

0.58 

0.88 

0.56 

1.21 

L2f 

0.58 

1.10 

0.51 

1.36 

0.63 

1.04 

0.41 

1.30 

Mean 

0.62 

0.97 

0.62 

1.34 

0.60 

0.91 

0.54 

1.24 


All units cm 
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Table B.5. ICESat SLR Residuals for 2011 Reprocessing POD 


Campaign 

Pass # 

Data # 

RMS 

Range 

Bias 

Radial 

RMS 
(data #) 
>60° 

Range 

Bias 

>60° 

Lla 

22 

330 

1.31 

1.01 

0.80 

n/a 

n/a 

L2a 

8 

272 

1.78 

1.03 

1.11 

1.82(10) 

0.53 

L2b 

3 

37 

2.34 

0.22 

1.26 

n/a 

n/a 

L2c 

23 

250 

1.57 

0.96 

1.07 

1.32(2) 

n/a 

L3a 

6 

48 

1.35 

0.78 

0.85 

n/a 

n/a 

L3b 

37 

999 

1.55 

0.97 

0.90 

1.50(13) 

0.76 

L3c 

68 

1860 

1.83 

0.90 

1.13 

1.52(51) 

1.16 

L3d 

191 

6278 

1.66 

1.06 

1.03 

1.56(142) 

1.43 

L3e 

85 

3146 

1.74 

1.20 

1.08 

1.60(65) 

1.82 

L3f 

216 

7741 

2.16 

1.15 

1.37 

1.75(311) 

1.52 

L3g 

241 

8803 

1.91 

1.08 

1.21 

1.88(346) 

1.52 

L3h 

222 

8561 

1.45 

0.89 

0.91 

1.45(304) 

1.38 

L3i 

235 

8270 

1.76 

1.01 

1.09 

1.42(267) 

1.09 

L3j 

115 

5103 

1.72 

1.05 

1.09 

1.57(225) 

1.37 

L3k 

99 

3813 

1.28 

1.01 

0.83 

1.48(127) 

1.44 

L2d 

86 

3308 

1.85 

1.21 

1.16 

2.09(113) 

1.71 

L2e 

178 

7113 

1.73 

1.16 

1.04 

1.58(279) 

1.40 

L2f 

53 

2290 

1.65 

1.19 

1.03 

1.94(54) 

1.66 

Total 

1888 

68222 

1.75 

1.04 

1.09 

1.58(2309) 

1.40 


All units cm 


Table B.6. ICESat SLR Residuals by Stations for 2011 Reprocessing POD 


Station 

Pass # 

Data # 

RMS 

Range 

Bias 

Radial 

RMS 
(data #) 
>60° 

Range 

Bias 

>60° 

Zimmerwald-I 

205 

5612 

2.11 

1.10 

1.33 

1.85(187) 

1.23 

Zimmerwald-V 

253 

8181 

1.73 

0.89 

1.08 

1.40(261) 

1.17 

McDonald 

86 

1147 

1.89 

1.36 

1.12 

n/a 

n/a 

McDonald2 

17 

316 

2.71 

1.50 

1.55 

n/a 

n/a 

Yarragadee 

555 

24847 

1.47 

0.98 

0.91 

1.64(697) 

1.42 

Greenbelt 

134 

5629 

1.76 

1.04 

1.01 

1.43(111) 

1.05 

Monument Peak 

144 

5072 

1.65 

1.00 

1.08 

1.48(213) 

0.98 

Haleakala 

27 

916 

1.67 

1.24 

1.07 

2.55( 6) 

0.00 

Graz 

227 

9186 

1.82 

0.94 

1.06 

1.36(257) 

0.99 

Herstmonceux 

171 

5729 

1.70 

0.98 

1.24 

1.61(523) 

1.30 

Arequipa 

23 

265 

1.48 

1.11 

1.02 

2.56( 9) 

2.68 

Hartebeesthoek 

46 

1322 

3.88 

0.80 

2.42 

1.84( 45) 

1.14 

Total 

1888 

68222 

1.75 

1.04 

1.09 

1.58(2309) 

1.40 


All units cm 
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